August 2021 


Foundations of bidomain models in electro cardiology 
Roberto Sudrez-Antola 


OMNIA Sciences/Technologies/Services 


Abstract: The equations of bidomain theory of electro cardiology are derived, from a suitable version of 
Maxwell equations, using an approach framed in the theory of averaged non-periodical fields. Micro, meso, 
and macro spatial scales are introduced, and a family of representative averaging volume elements is 
defined at the mesoscale. This volume averaging procedure is compatible with the histological and the 
anatomical details of the organ. It does not require the spatial periodicity assumed in multiple scale 
asymptotic expansions usually applied in homogenization problems. The macroscopic versions of charge 
balance, and Ohm law for the intracellular and the extracellular continuum are obtained by volume 
averaging. The effective conductivity tensors are derived solving a closure problem in a representative 
unitary element. A probabilistic approach framed in probability theorems of large numbers is used to derive 
the equation for membrane ionic current from the stochastic activity of channels at the microscopic level. 
An operational procedure suitable to define a sharp bidomain boundary from the fuzzy distribution of 
structural details and physical properties at the histological level are given. The bidomain equations are 
rewritten to introduce two new continua, both heterogeneous and anisotropic: a passive bioelectric medium 
and an active bioelectric medium. These media are coupled by two mechanisms: a surface coupling at the 
bidomain boundary and a volume coupling through a term that generalizes the activating function to 
anisotropic syncytia. A derivation of Maxwell equations for continuous media from an atomic scale version 
is briefly reviewed. A formulation suitable to describe the microscale processes in cardiac tissue is given 
and general boundary conditions at the interface between two media are stated. The compatibility of 
bidomain models with classical electromagnetic theory of continuous media is briefly discussed, 
considering two levels of averaging. 


(A) Introduction 


One of the most important features of the heart is its excitability. That implies its ability to respond 
both to internal (natural pacemakers) and external (artificial pacemakers) stimuli, producing and 
propagating the action potentials that initiate the heart’s rhythmic contractions during a whole 
lifespan. 


Bidomain models are mathematically convenient but highly simplified representations of 
excitable tissues, like myocardium [46]. Despite their limitations, these models have proven to be 
flexible enough to enable useful theoretical approaches to significant physiological, 
pathophysiological, and therapeutic problems [3]. 


To study the emergence and the propagation of action potentials in cardiac tissue, at least three 
main components of the electric current flow in the heart must be considered: 


-A flow in the shared intracellular space of cardiac myocytes. Although the plasma membrane 
encloses the whole cell, adjacent cell membranes of cardiac myocytes are electrically coupled by 
gap junctions. Due to this connection, ions and small molecules can flow from the interior of one 
cell to the interior of an adjacent one: cardiac muscle behaves as an electrical and chemical 
functional syncytium. 

The intracellular spaces of cardiac myocytes contain nuclei, mitochondria, sarcoplasmic 
reticulum, myofibrils, and other structured subcellular components, including membrane-less 
organelles like the centrosome, the P bodies, the nucleolus. All these components are embedded 
in the cytoskeleton and in a micro-trabecular network, which bound a watery solution of ions and 
small molecules. Figure 1 shows some intracellular structures in a ventricular myocyte and gap 
junctions related with intercalated disks. 
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Figure 1. Left: A sketch of some intracellular membrane bound organelles and myofibrils in a ventricular 
cardiomyocyte. Right: A sketch of three intercalated disks (on the right of the figure and below), an 
enlargement showing a set of gap junctions (right wing and above) and an even larger enlargement showing 
the connecting channels through the plasma membranes of two cells adjacent. 


-A flow in a shared extracellular space that is defined relative to cardiac myocytes. This 
extracellular space includes fibroblasts, blood capillaries and other components embedded in an 
extracellular matrix of collagen fibres. These fibres connect contiguous myocytes, and myocytes 
to contiguous capillaries. All these components bound a watery solution of ions and molecules. 
Figure 2 shows a sketch of some extracellular structures. 


Figure 2. Left: Interconnection of cardiomyocytes with an extracellular matrix of collagen fibres. Right: 
Blood capillaries adjacent to interconnected cardiomyocytes. 


-A flow through the excitable cell membranes of cardiac myocytes. These membranes separate 
the shared intracellular space from the shared extracellular space. The membranes present deep 
invaginations, the T-tubule. Through the intercalated disks, membranes establish mechanical 
(desmosomes and fascia adherents) and electrical connections (connexons in gap junctions) with 
the membranes of adjacent cardiac myocytes. 


A network of cardiac fibres results from the connections between cardiac myocytes. 

This network shows a high degree of regional specialization, structural and functional, closely 
related with the presence of different cell types: node cells, conduction cells and contraction cells 
[2] [21] [44]. These differences must be considered in the construction of mathematical models 
of the atria (a functional syncytium in themselves), the ventricles (another functional syncytium) 
and the propagation in the specialized atrio-ventricular connection. 

Figure 3 suggests the histological complexity of ventricular myocardium near the endocardium. 
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Figure 3 Partial histological section of ventricular muscle. Left side: longitudinal and transverse sections 
of cardiac muscle fibres formed by cardiac myocytes connected by intercalated discs. A venule and an 
arteriole can be seen in cross section. Right side: longitudinal and transverse sections of the conduction 
fibres of Purkinje. A cross section of a capillary can be seen, as well as subendocardial connective tissue 
and the endocardium. The sample was stained with Mallory-azan and is observed under a microscope at 
high magnification. Adapted from di Fiore’s Atlas of Histology with functional correlations [9]. 


In bidomain models, most of these morphological and functional details of cardiac tissue are 
ignored. Having in mind the above mentioned three main components of the electric current flow 
in the heart, the complexities of the actual tissue are replaced by three continua: intracellular, 
extracellular and membrane. 

Each continuum fills the whole spatial region occupied by cardiac tissue. The same spatial 
coordinate system is used to describe the three continua. At every point of the bidomain, the three 
continua are present, and the membrane continuum separates the intracellular from the 
extracellular continuum. 


This homogenization of the three geometric and physical elements that are necessary to describe 
the electric current flow in the heart (intracellular space, extracellular space, and membranes), 
allows us to construct a relatively simple mathematical framework. 

A system of coupled partial differential equations are posed in the bidomain. These equations 
describe an effective, macroscopic behaviour of the complex medium formed by cardiac tissue. 


The idea of homogenizing complex materials, introducing effective properties, can be traced back 
to the last 19 century. It was applied to flow in porous media, viscoelastic behaviour of materials, 
electric and viscous properties of dilute suspensions [6] [7] [50]. 

The origin of the effective medium laws and their effective parameters was mainly empirical. 
Later, in the 19 and mainly in the 20 centuries, these laws received several theoretical foundations, 
either deterministic or stochastic. These theoretical foundations are called upscaling. 

All deterministic upscaling methods have certain common features. 


They begin with microscopic equations my ;{a}) =( where 7 is a microscopic operator, VY 
is a microscopic field, and {a \ is a set of microscopic parameters. 


The upscaling operation (homogenization) (: . -) is applied to the microscopic equation. 


When the medium has microscopic scales widely separated from its macroscopic scales, the so- 
called localization approximation can be introduced: (aly; {ar} ))= M (yw ye B}) 


Misa homogenized (upscaled) operator, (y )is a macroscopic (homogenized) field, and { B } 


are effective (macroscopic) parameters. The resultant upscaled equation is M ((y ye B}) =0 
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The effective parameters can be determined experimentally. These parameters can be calculated 
also from the solution of a closure problem, posed in a representative region of the medium to 
decouple the upscaled behaviour from the microscopic details. 


Let now consider some historical milestones in the development of bidomain models applied in 
electrophysiology. 

In 1969, O. Schmitt proposed a description of certain physiological processes using the concept 
of interpenetrating domains [33]. 

In 1975, this idea was applied by Jack, Noble and Tsien to construct a three-dimensional cable 
model of biological tissues [19]. 

These authors described the electric current spread in a network of cells, with nearly the same cell 
size, and low resistance pathways connecting adjacent cells in the network. The resulting 
mathematical model corresponds to a homogeneous and isotropic functional syncytium, like 
myocardium or electrically coupled smooth muscle. 


However, the electrophysiological experiments done in the seventies, already showed that the 
electric conductivities in ventricular muscle, measured in the direction of the cardiac fibres, is 
different from the electric conductivities measured in directions transverse to the fibres, both in 
the intracellular and in the extracellular spaces [39]. 

Moreover, the anisotropy ratios (the quotients of longitudinal to transverse conductivities) in the 
intracellular spaces are generally different from the anisotropy ratios of the extracellular space: 
this is called unequal anisotropy [39]. 


In 1978, Tung [56] as well as Miller and Geselowitz [31], took all this into account and gave the 
standard mathematical formulation of bidomain models in electro cardiology. 

During the eighties, the nineties, and beyond, this formulation was significantly improved and 
applied to a wide spectrum of cardiological problems of physiological and medical importance. 
The bidomain equations were initially posed, and justified, based on intuitive and qualitative 
arguments. 

The proponents of the bidomain model did not attempt a derivation of the macroscopic equations 
from the known microstructure of cardiac muscle. 

They assumed that voltages and electric currents, as well as the effective parameters that appear 
in the bidomain equations, are macroscopic averages. 

Because, after averaging, detailed variations on a cellular and subcellular levels should be 
eliminated, a resolution capacity of a few cells was expected [38]. 


The first mathematical derivations of electro cardiology’s bidomain equations appeared in the 
nineties. 

Two different homogenization procedures were applied to filter out high frequency spatial and 
temporal variations, related with microscopic processes in cardiac muscle. 

In 1991, the method of volume averaging was applied by the present author to give a foundation 
to the bidomain models of electro cardiology [52]. 

In 1993, Neu and Krassowska applied an asymptotic two-scale method, to homogenize a periodic 
representation of the microstructure of a syncytial tissue [32]. 

In 1997, Krassowska applied this approach to a mathematical analysis of electroporation in 
myocardium [24]. ! 


' In 1998, the homogenization of the electric current flow in myocardium, done by the asymptotic two- 
scale method, already appeared in the first edition of the book by J. Keener and J. Sneyd [22], as a suitable 
foundation for a macroscopic approach to the propagation of cardiac action potentials. A detailed 
mathematical development of the asymptotic multi-scale homogenization method can be found in the book 


[5]. 
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The method of averaged fields presents certain advantages relative to the asymptotic multi-scale 
method, and it is intermediate between deterministic methods that require regular micro scale 
geometries and stochastic methods that operate with probabilistic ensembles to describe random 
microscale geometries [52]. 

Volume averaging procedure is compatible with the histological and the anatomical details of the 
organ. 

It does not require the spatial periodicity assumed in multiple scale asymptotic expansions usually 
applied in homogenization problems. 

Due to this versatility, in 1981, Nicholson and Phillips applied the method to model ion diffusion 
in the extracellular space of the rat cerebellum, drawing on the mathematical developments in the 
theory of transport processes in porous media [33]. This last work inspired the present author to 
apply the theory of averaged fields in cardiac electrophysiology. 


Our main objective here is to present a detailed derivation of the bidomain models of electro 
cardiology, using an approach framed in the theory of averaged fields. * 

We begin, in Section (B), with a brief review of standard bidomain models for electro cardiology. 
Then, in Section (C), three widely separated spatial scales are introduced in cardiac tissue 


(microscale or d scale; mesoscale or , scale ; and macroscale or L scale), jointly with a related 
set of time scales. 


The d scale is associated with the microscopic structural restrictions to processes in cardiac tissue, 


the /, scale is associated with the inevitable filtering processes of high spatial and temporal 


frequencies done by the measurement systems based in microelectrodes, and the L scale is 
associated with the measured values of the linear electrotonic parameters during an under- 
threshold polarization of myocardium. ? 

The associated time scales describe local dielectric and conductive relaxation processes at the 
microscale level and at the mesoscale level. The time scale associated with the macroscale level 
describe the time course of the foot of an action potential and related phenomena. 

In Section (D), a family of representative volume elements (REV) is defined at the mesoscale. A 
summary of the main results of the theory of averaged fields over the REV family is given. 

An error estimation of results obtained by the averaging procedure is given. The theoretical 
framework is applied to the volume average of the microscale fields and equations that describe 
the transport processes of a substance-like quantity (mass, charge, thermal energy, etc.). A 
macroscale version of a generalized Fick’s is obtained by volume averaging in a biphasic medium. 
Effective generalized conductivity tensors are derived from microscale processes, solving a 
closure problem by two different approaches. 


? In 2009 M. Perego presented, in the Politecnico di Milano, his Ph. D. thesis [34]. In this work, a concise 
derivation of the charge conservation in the intracellular and the extracellular continuum by volume 
averaging is given (the other equations of the bidomain models are postulated, as usual). Most of the 
research results obtained and discussed by the present author in 1991 [52] (even though fully available in 
Spanish, both in ResearchGate and in Academia), as well as some new results presented here, as far as the 
author knows, remained unpublished. 


3 In addition to these three scales, it is necessary to consider a fourth scale, atomic or a scale. In this last 
scale, the version of Maxwell's equations, in Physics, is called microscopic. Here we will say that it is an 


atomic scale version of those equations. Besides, the introduction of a fifth scale or 4, scale is necessary 


to make a volume average of the atomic level Maxwell equations and derive the equations of the 
electrodynamics of continuous media. 
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In Section (E) is briefly reviewed a derivation of Maxwell equations for continuous media from 
an atomic scale version. A formulation suitable to describe the microscale processes in cardiac 
tissue is given and general boundary conditions at the interface between two media are stated. 
The compatibility of bidomain models with classical electromagnetic theory of continuous media 


is briefly discussed, considering two levels of averaging. One level of averaging (4, scale) allows 
to pass from the atomic molecular version of Maxwell equations to the equations of the 
electrodynamics of continuous media. The other level of averaging (, scale) allows to obtain the 


bidomain equations as averages of the already averaged fields that appear in the electrodynamics 
of continuous media. 

The mesoscale averaged continuity equations and the mesoscale averaged versions of Ohm law 
for the intracellular and the extracellular continuum are obtained, in section (F), by volume 
averaging of the microscopic equations over the REV family. 

A representative unit element (RUE) of the tissue’s microstructure is introduced, and its 
relationship with a REV is considered. A procedure to derive the effective conductivity tensors 
solving a closure problem in a RUE is sketched. 

A probabilistic approach framed in probability theorems of large numbers is used to derive the 
equation for membrane ionic current from the stochastic activity of channels at the microscopic 
level. 

An operational procedure suitable to define a sharp bidomain boundary from the fuzzy 
distribution of structural details and physical properties at the histological level is given. 

In Section (G), the bidomain equations are rewritten to introduce two new continua, both 
heterogeneous and anisotropic: a passive bioelectric medium and an active bioelectric medium. 
These media are coupled by two mechanisms: a surface coupling at the bidomain boundary and 
a volume coupling through a term that generalizes the activating function to anisotropic syncytia. 
In section (H) some of the obtained results are discussed and a few conclusions and suggestions 
for further work are posed. 

In the Appendix of section (1), a dimensionless frequency domain version of Maxwell equations 
for electrodynamics of continuous media, formulated to be applicable to biological tissues, is 
analysed in terms of a regular perturbation method based in small dimensionless parameters that 
reflects tissue’s properties. 

Finally, in section (J) a list of main references is given. 


Let us begin summarizing the standard formulation of bidomain models. 
Before doing this, it seems advisable: 


(a) Introduce a convention in relation to numbering and the reference to equations in the same 
section of the work or in another of the sections. 

The equations will be numbered in each section. So, inside the same section an equation will be 
referred to only by its number, but in another section, it will be referred to with the letter 
corresponding to the section to which the equation belongs, followed by the equation number 
there. 


(b) Describe the bibliographic convention used in this work: some of the references are given in 
full, forming part of footnotes, while others, those considered of greater importance, are presented 
at the end of the text, in the section of bibliographic references. There they appear in alphabetical 
order, preceded by a number in parentheses: []. 


(c) Introduce a convention on the mathematical symbols used and a short glossary of some of the 
main symbols used. 
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The mathematical symbols of set theory, calculus, vector, and tensor (dyadic) analysis are the 
most often used. 4 

A possible exception is the vector product and the curl: instead of < we use the symbol A. 
Section (E) begins considering condensed matter with a level of resolution at the scale of 
Angstroms (a scale). 

At this level of resolution, both electrons and nuclei can be thought as point charges. 

Then, Maxwell-Lorentz equations in Euclidean space are taken as the starting point to summarize 
the averaging process that ends in the equations of electrodynamics of continuous media. 

Due to atomic movements, at the a scale the boundary of any system of condensed matter is 
blurred. Because of this fuzziness, the different instantaneous distributions of the point charges 


in the whole Euclidean space R° is considered. 
The position vectors of the points of the condensed matter system, when considered at the a 


scale, are represented by fF € R’. 
Now, let us make an average at the d , scale (a hundred Angstroms) as will be sketched in section 


(E). Instead of a swarm of particles, we obtain a continuous medium filling a whole region. 


The properties of this medium change significantly only at the d scale (with a level of resolution 
at the scale of WM), 


When referring to the region occupied by the cardiac tissue described in the d scale, the volume 


of space occupied by the biological tissue is represented by a region QC R? , assumed to be well 
defined and with a regular boundary surface. 


The intracellular space is represented by Q ;, the extracellular (interstitial) space is represented 
by Q, -, and the membranes separating the intracellular from the extracellular space is represented 
by the surface A,,. Then QD=Q,UA,UQ, . 

The extra tissular space, when considered in the d scale, is represented by 2. The interface 
between the cardiac tissue {1 and the extra tissular space (is represented by OQ. The position 


vectors of the points of the cardiac tissue or of the extra-tissue region, when considered at the d 
scale, are represented by y (so yeQ ; is a point belonging to the intracellular space, y € A, is 
a point belonging to the cell membranes, etc.) 

However, if they are strictly considered as points in three-dimensional physical space, r and y 
could be the same point, represented with two different notations, depending on whether the 
averaging at the scale @, is already done (y ), or it is not ce ). 

When referring to the region occupied by the three overlapping continua (intracellular, 


extracellular and membranes), obtained after averaging in cardiac tissue in the (, scale, it is 
represented by the symbol ‘t. The extra tissular space, when considered after averaging in the ty 


scale, is represented by ‘8, and the interface between these two regions, after averaging in the i 
scale, is represented by OR. 
The position vectors of the points of the cardiac tissue or of the extra-tissue region, when 


considered after averaging at the (, scale, are represented by x. 


+ A useful compendium of mathematical and cross-disciplinary terms used in the sciences and engineering 
can be found in Previato E (ed) (2003) Dictionary of applied math for engineers and scientists, CRC Press, 
boca Raton, Fl USA. 
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Now, due to the construction of the bidomain models, each point X belongs simultaneously to 
three continua: intracellular, extracellular (interstitial), and membranes. 

However, if they are strictly considered as points in three-dimensional physical space, y and X 
could be the same point, represented with two different notations, depending on whether the 


averaging at the scale (, is done, or it is not. 


Each point X is the centroid of a representative elementary volume (REV) B (x) , while B , (x) 
is the subset of B (x) that belongs to phase /. The indicator function of phase // in the REV is 
I =| eae 

v 0 if ye¢B, (x) 


A scalar field w (t, y ) defined in Q or ina part of Q is called a local point field. 


defined to be: 


The extrinsic phase average of the field in the point X € R is defined relative to the whole REV 
volume, being y, (¢, y) =y (¢, y) I, (y) and V (B(x) the volume of the REV: 


a Se ae thes 
(v,)(t%) “TBs (7.9) 1,(y) av (9) 
The intrinsic phase average of the field in the point X is defined this way: 
Pr hee a 1 4 4 4 
(Yu) "TB ®) B(z) (1, ¥) 1, (¥)av(¥) 


Phase averages of a given field, either extrinsic or intrinsic, will be called local averages of the 
fields. 

Orthogonal cartesian coordinates will be used in all cases. 

Everything said about the averages of scalar fields can and will be applied to the cartesian 
components of vector and tensor fields. 


(B) Asummary of standard formulation of bidomain models of electro cardiology. The 
monodomain model of the membrane continuum. 


Three superposed continua are supposed to fill the whole spatial region occupied by 
myocardium: the intracellular, the extracellular and the excitable membranes. 


The electric potential in each continuum is represented by @(t,x): intracellular @ (t,x), 
extracellular , (t,x) and membrane ©,, (x): Here ¢ is the time instant and 


X=xE +y a a4 é, is the position vector in orthogonal cartesian coordinates with unit vectors 


= 


é.,€,,€.. At each point and at each time instant it is assumed that the intracellular electric 


xo Vvyr&z 


current density J, is given by 


I= «(Zoe «(22 }e-s.{ 2 Je (1 a) 


=> 


The extracellular electric current density J 


is assumed to be given by 


e 
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_ a 2 a a 7 
A se (Ze, lee, (Zo, }e -s.{2 Ay 0, Je (1b) 


By a suitable selection of coordinate axes, the electric conductivities §;, , intracellular g;, and 
extracellular &,, correspond to fibres direction, while 8;,, 8.,, 8),, 8-, are transverse 
conductivities. These conductivities may be functions of the position in the bidomain. 


Also, due to local changes in blood volume, interstitial space fluids or intracellular medium, the 
conductivities may be functions of time. 


The conductivities g; and g, are assumed to be related with the volume fractions of the 
intracellular domain f; , and of the extracellular domain f,, and with their intrinsic 


conductivities O; and O,: 
Bi, = 5, x 8 =F, 5, Siz =f, 0;, Sex =F. Cex Bi = f0,5 802 =F,%e; (2) 


However, electric current densities J (t,x) and J : (t,x) are defined as the electric charge 


crossing a unit area in a unit time, and not referred to volume. As consequence, equations (2) 
should be formulated in terms of the area fractions taken on three mutually orthogonal cross 
sections in each point of the bidomain. 


If Tek is the area fraction of intracellular space in a cross section orthogonal to the x axis, and 


so on, we would have these relationships, instead of the ones given in (2): 


Six = five ix Siy = Te, Siz =e ey 07. (3a) 


Sex = Fe yeOes Sey =f, xc ey Sez =f xyFer (3b) 


Besides unequal anisotropy, there is another fact that is observed in ventricular myocardium, 
related with its complex three-dimensional spatial organization. The elongated cells are 
connected, mainly end to end, forming a hierarchical structure of sheets (lamina) of cardiac fibres. 
The sheets, separated by layers of connective tissue, rotate smoothly from the endocardium to the 
epicardium (Figure 4). 


Figure 4. Left side: a planar sketch of a ventricular cardiac myocyte with its connections through 
intercalated discs to adjacent cells. Right side: rotation of muscular sheets from the endocardium to the 
epicardium. 

To take this geometry into consideration, an orthogonal curvilinear coordinate system can be 


introduced, with unit vectors &,() parallel to the local fibre direction, &,(X) tangent to the 
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lamina there, and é,(x) orthogonal to the lamina. Then, the generalized Ohm laws given by the 


equations (1) and (2) can be rewritten in curvilinear coordinates. 


The potential in the membrane continuum ©, (t, x) is assumed to be given by: 

®,,(t,¥)=®, (t,x)—®, (t,x) (4) 
If there are not sources of external current, neither in the intracellular nor in the extracellular 
continuum, and if there is not charge accumulation in these continua, the current leaving 


(entering) the intracellular domain must be equal to the current leaving (entering) the extracellular 
domain. 


Then the divergences (represented by V-(...)) of the intracellular and extracellular current 
density fields verify: 

V-JF,=V- I =7_J, (5a) 
Here 7,, J, is the current per unit volume that crosses the membrane continuum (positive when 
directed from the intracellular towards the extracellular domain). 
The parameter 7,, gives the local membrane area per unit volume, while J,, gives the local 


electric current density per unit area of membrane. > 
From equation (5 a) we obtain the following quasistatic charge conservation equation in the 
bidomain: 


V-(J,+4,)=0 (5b) 


In each case, to be able to detail the components of the local electric current density per unit area 
of membrane, a mathematical model of the membrane must be defined. 


All of them assume that J,, is a function J,,(® W) of the field of membrane voltage 


m? 


OQ, (t, x) and a n-dimensional vector field W (t, x) defined in the membrane continuum. 


The scalar components of this vector field represent the activation and recovery variables. 
The nature and number of these scalar components, as well as their interrelation with membrane 


voltage © , are defined in each membrane model. 


In a detailed analysis, different models must be applied to Purkinje fibres, sinoatrial node cells, 
atrial muscle cells, atrioventricular node cells, and ventricular muscle cells. © 


5 From equation (4) we obtain the following quasistatic charge conservation equation in the bidomain: 
V-(J, +J,) -0 When there are sources of external current, either in the intracellular or in the 
extracellular continuum, this equation must be modified as shown below. If the electric flow in the tissue 


is produced by external electrodes only, and the if the membrane currents in the sarcolemma of the cardiac 
myocytes can be neglected, at least as a coarse approximation, then the electric current flows through the 


extracellular space in quasistatic conditions: V-J Fim 0 


° A summary of cardiac muscle cells models can be found in a chapter Membrane models, written by A. 
Varghese [58]. In the already mentioned books [5] and mainly in [62] it can be found a detailed exposition 
of membrane phenomena and membrane models. 


10 
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The function J »(®,, W) is also defined in each membrane model, although it always involves 


an ionic current and a displacement current component. 


m? 


If the electric consequences of membrane invaginations can be ignored, then J »(® W) is the 


sum of a density of displacement current and a density of ionic currents: 


J,( pW) =Cy <u ( Ppa H?) (6) 


C. is the electric capacity of the unit area of membrane, and the density of ionic current is 


represented by J,,,,,;. ( Q,. W ) 


If the consequences of membrane invaginations must be considered, as well as the effect of 
internal membranes and the narrow clefts between adjacent cell membranes, the simplest 
mathematical model that allows the description of the membranes is a four components circuit. 
Such a model could be a nonlinear extension of the mathematical model proposed by Falk and 
Fatt to interpret the linear electrical properties of skeletal muscle fibres observed by intracellular 
microelectrodes [13] [19]. ’ 

However, although it can be of interest in a study of the electric behaviour of ventricular muscle, 
we don’t consider this and other extended mathematical models in the present research. 


The kinetic of the activation and recovery variables in J »(®,,,.W) is described by a nonlinear 


system of ordinary differential equations: 


Ouse & 2 
a W = F(nW) (7) 
Besides the bidomain that represents cardiac tissue, a third dominion is introduced: a volume 
conductor surrounding it. In general, the electric current flow in this volume conductor is 
described by a macroscopic version of Ohm’s law, in terms of a single macroscopic scalar 


potential ©, and the corresponding conductivities. To simplify, sometimes the external volume 


conductor is assumed to be homogeneous and isotropic. 


At the boundary between the cardiac bidomain and the surrounding continuum, the continuity of 


the extracellular potential ©, with the external potential ) is postulated. 


The flow of current from the intracellular domain to the external domain, through the membranes 
of cardiac cells located at the boundary between the cardiac bidomain and the surrounding 
continuum, is assumed to be negligible. 


7 Tn a four element model a suitably defined membrane current density is in parallel with a current density 
through an in-series combination of a linear resistance R , of a unit area with a capacitance C , per unit 


area. If c¢, =R,C, 1s a relaxation time, then the global density of membrane current J, ‘ resulting from 
all the effects due to the membranes is related with membrane voltage through the equation 


é oe 
tng tig OG Pa (Cc 


S m 


+C, s ©, +7,—J LJ A linearized version follows 
t 


S Ct tonic, g : ionic, g 


from the nonlinear equation by the substitution 7 = = “(0 -®,) , where ®, is the rest voltage and R 
tonic ,g m m 
m 


is the resistance of a unit area of membrane [10] [11] [13]. 


11 
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As consequence, all the external current must enter or leave the cardiac bidomain through the 
extracellular continuum. 


Nevertheless, the possibility of current injection with a microelectrode located in the intracellular 
or in the extracellular spaces cannot be discarded. 
jeiae 9 Qfepresents the current density per unit volume injected by a microelectrode in the 
extracellular space, then: 
V-I,=7,5 41, (8) 
From (5 a) and (8): 
V-(S,+5,) =I (9) 
When the microelectrode is located at the point X, of the cardiac bidomain, then: 
1,(t.8) =i(t)5(#-%,) (10) 
The electric current injected is represented by i(t ) and 6 (x- X,) is a three-dimensional Dirac’s 


delta function (in fact, a distribution). 


Let us now reformulate the equations of the bidomain models of myocardium, using a tensor 
formulation, free of coordinates. For that purpose, we introduce conductivity tensors, G, for the 
intracellular domain, and G, for the extracellular domain. 
The generalized tensorial versions of Ohm’s laws are: 
J,=-G,V®, (11) 
J,=-G,V®, (12) 
From (4) it follows @ i= oO, + O, and from this last equality and from (5) and (9) the following 


coupled bidomain equations: 


InJn =V (GVO, )+V-(G-V®, } VEER (13) 


V-(G:Vo,,)+V-((G,+6, )-V,)=1, VEER (14) 


When the electric consequences of intracellular membrane invaginations and the possible 
ephaptic coupling between adjacent membranes of different cardiac myocytes can be neglected, 
then to equations (12) and (13) we can add equations (6) and (7) for the electric current density 


through the membrane continuum J, (,,.W). 


To complete the bidomain model, both initial and boundary conditions must be introduced for the 
fields D,, and ®,, and the initial values of the vector W of channels activation and recovery 
variables. 


Let us represent the region occupied by the cardiac bidomain by}, its boundary by OR and the 


region occupied by the external volume conductor by ¥ ,. 
Introducing a tensor conductivity G, for this last conductor, being @, the electric potential, and 


J, the electric current density there, the generalized Ohm’s law in this external continuum is: 


J, =-G,-V®, (15) 
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We assume that either ®,(t,x) or J,(t,x) are known in ). 


The unit normal to each point of the boundary, directed outwards from the region 8 is represented 


by 1. Then, if the flow of current from the intracellular domain to the external domain, through 
the membranes of cardiac cells located at the boundary, is assumed to be negligible, we have: 


iJ, =0 So, from 2, =®,, +, and from (11): 
i-G,V®,=-i'G,V®,  VXEOR (16) 


As consequence, all the external current must enter or leave the cardiac bidomain through the 
extracellular continuum: 


iJ,=i-J, 2 wG,VO,=-H-J, VXEOR (17) 


When G, =k G, , being k a dimensionless constant, the anisotropy is called equal. Introducing 


ay are 
an effective electric conductivity G = “cer G, , equations (13) and (14) reduce to one equation: 
+ 
In Jn =V (GVO, | VEER (18) 


When the cardiac bidomain is electrically isolated through is boundary, from (17) we obtain 


n- G, -V®, =0 and as consequence, from (16) it follows that: 


i-G,V®, =0 VxEOR (19) 


The equation (18) jointly with the boundary condition (19) are known as the monodomain 
approximation, formulated in terms of a membrane continuum. 


Four main problems can be distinguished in relation to bidomain models: 


(1)-To justify the introduction of the three superposed continua, including the fields that appear 
in the equations of bidomain models, considering the complexities of tissue’s microstructure at 
the histological level. 

(2)-To obtain, unless a certain error to be estimated, the bidomain equations 


ayy ,;=V: J -=X,,/7,, from the equations formulated at the microscopic scale. 


This should include a derivation of J »(®,,,.W) and the equations (6) and (7), considering 


microscopic details of the excitable membranes. 


(3)-To derive J ,= -G, ‘V®, and J 2= -G, ‘V®, from the equations of electromagnetism and 


transport processes in the microscale unless a certain error to be estimated. 

The effects due to irregularities of the space-time distribution of the fields and the complexities 
of boundary conditions originating in tissue’s microstructure, should be considered. 

This derivation should include a justification of equations (3a) showing that, unless a certain error 
to be estimated: 


(4)-To give a constructive definition of the interface between the bidomain and the extra tissular 
domain. 
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To apply the (upscaling) method of volume average in a derivation of the bidomain equations, 
giving a solution to the problems summarized above, as we said before, we must introduce three 
spatial scales (microscale, mesoscale, and macroscale) in cardiac tissue. 


Besides, the compatibility of bidomain models with classical electromagnetic theory of 
continuous media should be discussed, considering two levels of averaging. 


One level of averaging (4, scale) allows to pass from the atomic version of Maxwell equations to 
the equations of the electrodynamics of continuous media. 
The other level of averaging (£, scale) allows to obtain the bidomain equations as averages of the 


already averaged fields that appear in the electrodynamics of continuous media. 


We need five spatial scales: a (to formulate the atomic version of Maxwell equations), 4, (to 


obtain the equations of the electrodynamics of continuous media), d (to describe the structure of 
cardiac tissue at the histological level, and to describe the spatial variations of the electromagnetic 


field and the electric and magnetic properties of the tissue at this level), (, (to apply the theory of 


averaged fields in order to construct the bidomain model), and L (the scale of variation of the 
fields and tissue properties that appear in the bidomain models of electro cardiology). 

So, let us consider the problems that appear in relation with the definition of these scales, 
beginning with the last three of them. 


(C) Measurements, scales, and systems 


The presence of highly resistive cell membranes (and highly resistive membrane less structures 
in the intracellular and in the extracellular spaces) introduces microscopic constraints to the low 
frequency (less than 10 kHz) components of the electric current flow in the much less resistive 
intracellular and extracellular spaces. 

At high frequency (higher than 100 kHz), the impedance of the membrane is near zero (due to the 
capacity in parallel to the ionic channels), and the membranes behaves as short-circuited. 


However, bidomain models are constructed to describe low frequency bioelectric phenomena. 
So, the membrane constraint arguments can be applied in this case. 


Membranes are natural boundaries, both from the electric and the morphological points of view: 
it is possible to take the membranes as phase boundaries between an intracellular phase and an 
extracellular phase. 

The left side of Figure 5 shows an example. Phase a can be considered as a simplified 
representation of the intracellular space and phase f° can be seen as a simplified representation 


of the extracellular space. 
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Figure 5. Left side: Sketch of a two-phase medium. Right side: Definition of the local diameter in any point 
P interior to a phase, either a or £ . An open ball is constructed. The ball is included in the phase to which 


P belongs, has the point P as one of its elements, has the maximum possible radius 5( P) , and centre in 


C(P)- 


A local diameter can be defined at each point of these two phasesa and /, following the 


procedure proposed by A. Scheidegger to define local pore diameter in a porous medium [50]. 
Let us consider a not too small region of tissue (with linear dimension of the numerical order of 
1 mm). 

For each point P belonging to the interior of a given phase, a spherical open neighbourhood 
ELC 6 ) of centre C and radius 6 is identified, such that B(C;6 ) is the open sphere of 


maximum radius included in the phase to which P belongs. 


This is shown in the right side Figure 5 above, for a selection of four points F and P, in phase 
B, P, and EF, in phase a . The centres and radii of the balls are C, and O, for the chosen point 


Ey , being k =1,2,3,4. 
Then we take the ball’s diameter as the local diameter of the phase at point P. 
The mean value of these local diameters can be selected as an estimation of the microscopic scale 


of a given phase. ° 


The maximum of the resulting two mean values (intracellular d ; and extracellular d, ) is taken 


as an estimation of the macroscopic scale d for the whole tissue. 

According to this definition and considering the geometric counterparts of the physical 
restrictions due to the membranes, d follows from quantitative morphological details of the 
tissue. ° 


8 The reason for not taking the maximum of the local diameters corresponding to one of the two phases, 
instead of their mean value, is that a very large local diameter can appear that behaves like an outlier relative 
to the radii distribution. 


° An in-silico reconstruction of the microscopic geometry of a tissue, like cardiac tissue, is now possible 
applying the methods of quantitative morphology [3] [5] [29] [57] [59]. This topic is briefly considered in 
section (D). 
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An associated time scale to the microscale d, is the polarization relaxation time in the intracellular 


space. 
To obtain a coarse estimate of this time scale, we begin with a simple 1-D model of a cylindrical 


volume of cross section A and length d,, composed by a in series combination of an area A of 


purely capacitive membrane with electric capacitance per unit area C,,, , and an ohmic resistance 


R= om with intracellular representative conductivity O, . The polarization relaxation time in 
on 


L 


this case is: 


r=—"d, (1) 


A polarization relaxation time 7, is obtained for the extracellular space by substituting in (1) the 
representative conductivity for the intracellular space O, by a representative conductivity of the 
extracellular space O,, and the microscale d, by the microscale d, . 

The other two spatial scales, mesoscopic and macroscopic, are closely related with 
electrophysiological measurements done by microelectrodes. 


As an example of a possible electric measurement system employed in vitro in electrophysiology, 
let us consider a microelectrode inserted inside a cell, as suggested in Figure 6. 


— 


Figure 6. This sketch shows a glass microelectrode inserted in an intracellular space and connected to an 
amplifier A. The amplifier is connected to a voltage recording unit Ey. An auxiliary electrode, grounded, 
is shown in the physiologic solution that surrounds the measured cell. 


The microelectrode is connected to one of the terminals of an amplifier, while an auxiliary 
electrode, grounded, is connected to the other terminal. The output of the amplifier is connected 
to a voltage recording unit. 


The voltage signal y(t ,x) registered by the instrument can be considered as the result of 
progressive transformations undergone by a hypothetical input signal y(t ae in principle a 


field varying in space and in time. Here X is the position of the microelectrode inside the cell, 
assuming it can be represented by a point in space. 

These transformations depend on the characteristics of the electrodes, the amplifier, the recording 
device, and the connection elements between the subsystems of the measurement system 
(including possible effects related to electromagnetic compatibility). 
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For the purposes of the present discussion, and neglecting for the moment, random noises, the 
link between the assumed input and the measured output can be represented by a linear, causal, 


stationary and continuous operator with kernel K (t, x) defined in a small space region B (x) 


centred at X , and in the time interval [t —At, t] 70 


v.(t.3)= f atl [yg K (0.8). )aV (2) 2) 


t-At 
The input signal, as such, is not directly accessible through a measurement process, so the 


postulated K (t ; x) cannot be experimentally determined. 


However, we can conceive a sequence of measurement systems, such that: 
a) In each one of them Ar is lower than the value of Ar in the preceding system and greater than 


At in the following system. 
b) In each of them the region B (x ) is included in the corresponding region of the preceding 


system and includes the corresponding region of the following system, the dimensions of these 
regions decreasing when passing from one measurement system to the next in the sequence. 

c) The confidence interval, for a given level of statistical significance, corresponding to repeated 
measurements of the same variable carried out in one of the systems, always includes the 
confidence interval, with the same level of significance, obtained for the same variable, measured 
in the following system of the sequence. '! 


Then we can relate the field measured by a given system to the field measured by a more precise 
and accurate system, taking this last signal as a reference. 

In any case, the output signal of a measurement system can be modelled as a weighted space-time 
average of an input signal. 

In practice, this averaging includes the damage region that always appears when the 
microelectrode is introduced in the intracellular or in the extracellular spaces of cardiac tissue. 
And if the electrode is a glass microelectrode with a tip of a few tenths of a micrometre in 
diameter, the interaction with the biological tissue takes place through a small but finite liquid 
junction region. 

All electrodes measure voltages by averaging over a characteristic mesoscopic interaction region 
and filtering out details that appear at the microscopic scale. 

The mesoscopic spatial scale @, must be chosen at least as extensive as the region of interaction 


of the electrode with the cardiomyocyte’s interior or with the extracellular space. 

Besides this local mesoscopic interaction region, another consequence of the interactions between 
electrode and tissue is the appearance of a lower limit for the distance between points at which a 
variable such as voltage can be measured, for example, within the same cell. 

In addition to the inevitable error in voltage measurements, the micromanipulator that allows the 
microelectrode to be located within the tissue, always presents an error in position that must be 
considered. 


10 A random noise could be added to the output signal in equation (2). 


'! This is known as external convergence of measurements. External convergence allows us to define a 
physical quantity by referring to a class of measurement methods, rather than just one. This justifies the 
current use of the term physical quantity without making explicit the different methods that could be used 
to measure it [30]. 
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Under these conditions, the electrophysiologists can measure mesoscopic averaged voltages at 
different points of, for example, a steady polarized intracellular or extracellular space in a sample 
of myocardium. The source of polarization could be an injected current through a point electrode. 
As the voltage pattern remains well below the threshold of the action potentials, the polarization 
is electrotonic. 

The space constants appear as parameters in 1-D, 2-D, or 3-D mathematical models [28]. 

The parameters of a model that describes the specific experimental situation, are adjusted using 
the experimental results. Thus, the space constants of the sample of heart tissue are obtained. 

All other conditions being equal, the space constants corresponding to these tissue currents 
increase when the number of dimensions (in which the electrical currents are distributed in the 
tissue) decreases from 3 to 1. Results between hundreds of micrometres and a few millimetres are 
obtained [28] [38] [39]. 


The macroscopic spatial scale L of the cardiac tissue is defined as a distance through which the 
locally averaged voltage profiles and other related fields undergo significant variations. '” 
The space constants 4, can be used as bidomain macroscopic spatial scales. 


A coarse estimation can be obtained working with a mathematical model of a straight, cylindrical 


cardiac muscle fibre, of diameter d , uniform internal conductivity o, , and a membrane polarized 


near the rest voltage ©, > With resistance R.,, of the unit area of membrane and capacitance 


C, per unit area. The following cable equation is obtained under these restrictions: ' 


re) sO: 2 
loi nH Ou Pea) = Nes nae +A, 55 ®, (3) 


In equation (3) the space constant (macroscopic length scale) is: 
1 
Ay = [Rod (4) 


The time scale associated with the macroscopic length scale is the membrane’s time constant T,, 
FRC. (5) 

Let us summarize what was previously discussed in an example of the use of three spatial scales 

in the interpretation of the results of an electrophysiological experiment carried out with 

microelectrode. 

In a thought experiment three voltages are measured, in steady state, positioning a microelectrode 

in the intracellular space of three adjacent ventricular myocytes, members of the same myocardial 


"2 The definition of scale of a function and the concept of orthodoxy of function’s behaviour is carefully 
developed in the classic applied mathematics book by Lin and Segel [27]. 


'3 The electric current flow in the intrafibrillar space is given by: d; (t,x) = a a? Joy = o. (t, x) 
x 


The electric current flow through the unit length of fibre membrane is given by [ ‘i (t xX ) =nmdJ 24 (t xX ) 


The electric current per unit area of membrane polarized near the rest voltage can be described by the linear 


formula Jin = Crs ®, a2 1 (®,, ©, i) (positive outwards from the intrafibrillar space). 
t : 


m 


Besides [, (t,x) = -<1, (t,x) and , (t, x) =, (t,x)+®, (54) being ©, (t,x) the extracellular 
x 


voltage in the neighbourhood of the fibre membrane. From these relationships and the definitions (4) and 
(5) it follows equation (3). 
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fibre. The distance between these electrodes can be taken to be 100 zm, considering the length 


of ventricular myocytes. 
The position-voltage pairs obtained after processing the data appear as filled circles in the three 
following figures. 


Figure 7 suggests how the electric potential field ; vary as a function of the distance x. 


This field pattern can be thought as obtained, measuring simultaneously in every point of the 
sample, with a measuring system with greater field resolution capacity than a system one of whose 
transducers is a microelectrode. 


&, (x) 


0 


x 


Figure 7. Intracellular electric potential field 0, as a function of the distance along the fibre. 


Some noise should be expected in the signal, but to see it we should have measured during an 
interval of time. 

From now on we assume that the noise is filtered out. As consequence, the small-scale variations 
observed in the curve of Figure 7 can be due to small scale inhomogeneities in the electric 
properties of the intracellular space. 


The two steps of voltage drop correspond to the gap junctions present in the two intercalated disks 
that connect adjacent cells. 

The three sections of the voltage vs. position curve, where the average slope of the voltage is 
lower than the slope in the intercalated disks, correspond to the intracellular spaces of the three 
cells considered. 


Figure 8 suggests how the electric potential field ®, , measured by a glass microelectrode would 


i,a 
vary as a function of the distance x . Now, the measuring process implies a spatial averaging in 
a region around the electrode tip (damaged volume and junction potential) 


BS, y lu’ 


0 = 


Figure 8. Electric potential field G(x ) measured by a glass microelectrode. 


19 


August 2021 


The high frequency spatial components of the curve that appears in Figure 7 are filtered out. 
However, the two steps of voltage drop correspond to the gap junctions present in the two 
intercalated disks that connect adjacent cells, still appear in the voltage pattern. 


Figure 9 sketches the voltage (®,) averaged in the mesoscopic scale, as a function of distance. 


€B.>/x) 


Oo = 


Figure 9. Electric potential field (®,) averaged in the mesoscopic scale, as a function of distance x . 


Mesoscopic averaging eliminates spatial frequencies that averaging carried out with the 
microelectrode does not. In this example the voltage drops associated with the gap junctions 
disappeared, but this is not always so. 

Figures 7, 8 and 9 show that the three measurements done in the steady state thought experiment, 
taken in isolation, are compatible with each one of the three assumed curves. 


If the electrophysiological context is not a steady one, a dynamic analysis must be done. 

It is possible to simplify the dynamic analysis of a dissipative system like myocardium, with 
widely separated time scales, applying the following two principles to link time scales of different 
orders of magnitude. 


A reference time scale must be selected first. Then: 


1-The variables belonging to processes with time scales at least an order of magnitude greater 
than the reference time scale, can be considered as frozen. 


2- The variables belonging to processes that evolve with time scales at least an order of magnitude 
smaller than the variables of reference, after a short transient, can be considered as relaxed to 


equilibrium with these variables, evolving in the so-called outer time scale. 


For a process at the bidomain level, the time scale can be estimated applying formula (5). 


With R,=10kQxcm* and C, =1 ase we obtain: 
cm 


T,, =10ms 


The polarization relaxation time in a conductive intracellular space bounded by purely capacitive 
membranes, can be estimated by formula (1). For a representative intracellular conductivity 


mS 
o = 4— , and assuming, only for the moment, that d, is equal to a possible ventricular cardiac 
cm 


, then formula (1) gives: 


LF 
=1-. 
cm 


myocyte diameter, that is 10 um, and C,, 
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7, =0.25 ps 
So, the dimensionless number that results from the quotient of the two considered time scales is: 
tT, 0.25 us r 
Fi = 95x 10° 
tT,  LOms 


Even if we suppose that d, is equal to a possible ventricular cardiac myocyte length, that is 
100 um, a very small estimation of the relaxation time (2.5 szs ) is obtained. Similar results are 


obtained for the polarization relaxation time in the extracellular space. 


This example suggests that we could apply the second principle to link the time scales of 
microscopic processes with the time scales of macroscopic processes. 

According to this principle, when the last time scales are taken as reference scales, the polarization 
variables of the tissue are permanently in equilibrium with the variables that evolve at the level 
of the bidomain (macroscopic) time scale. 


From equation (4) we obtain the following estimation of the macroscopic length scale of the 
bidomain: 
A,, =|lmm 


Now, let us discuss the construction of a suitable mesoscopic scale for volume averaging of fields. 


(D) Mesoscopic scale, volume and area averaging, generalized macroscopic Fick’s law, 
closure, and other related topics 


Here, the methods of volume averaging, applied to transport processes in a model of biphasic 
medium (with phases @ and £ ) is described. These processes are not restricted to electric charge 


transport and the model of biphasic medium is not restricted to a model of cardiac tissue. 

A family of representative volume elements (REV) is defined at the mesoscopic scale in the 
biphasic medium. A summary of the main results of the theory of averaged fields over the REV 
family is given. An error estimation of results obtained by the averaging procedure is suggested. 
Then, the theoretical framework is applied to the volume and area average of the microscopic 
fields and equations that describe the transport processes of a substance-like quantity (mass, 
electric charge, thermal energy, etc.). A macroscopic version of a generalized Fick’s law is 
obtained by volume averaging in the biphasic medium. Effective generalized conductivity tensors 
are derived from microscale, generalized conductivity tensors, solving a closure problem by two 
different methods. 

In the theory of volume averaged mass, momentum and heat transport processes in porous media, 
the mesoscopic scale is often chosen as the geometric mean of the microscopic d and the 


¢,=JLd (1) 


When the microscopic scale is enough separated from the macroscopic scale: 


li 
Ke ceed oe (2) 
CONES E 


macroscopic L scales [33] [60]: 


If the restrictions imposed in (2) are obeyed, then simultaneously d « ¢ , and t .< Lx, 


To get closer to this and other possible definitions of the mesoscopic length scale, for each point 


P of a biphasic medium, we define a cubic volume region B (Gz ) of length’s side es 
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This is called a representative elementary volume (REV). The point P is the centroid of this 
volume element. 

The anisotropy of the cubic volume favours some spatial directions, which does not occur in the 
case of a spherical volume. Besides, both give a steplike behaviour for the averaging process 
when crossing the boundary from the interior to the exterior of the averaging volume, because the 
averaging procedure gives the same weight to all the values of the local point field, with 
independence of the position relative to the boundary of the REV. 


There are other, more sophisticated definitions of volume average in complex media [7], but we 
will not consider them in the present research. However, see next Section (E) Quasi-steady 
microscopic electric flow in cardiac tissue. 

Furthermore, we will introduce this restriction to the possible cubic REVs: 


(R1) The size (given by the mesoscopic length scale Ly shape and orientation of the 
representative elementary volume is independent of the position of the centroid in the medium. ' 


So, the volume of the cube must be the same for all the REVs: V (B (P)) = (3 


In fact, instead of a cube, the REVs could be a set spheres or any other set of regions (congruent 
by translation without rotation) suitable to apply the averaging procedure. 


Besides, to obtain useful averages it is necessary to introduce two additional restrictions to the 
mesoscopic length scale @, [51]: 


(R2) d<¢/,  Assures the filtering out of microscopic fluctuations in the fields and allows the 


derivation of approximate macroscopic equations by averaging the microscopic equations. 
(R3) ¢,<ZL_ Assures a coherence condition: the local average of the local average of a field 


can be equated to the local average of that field. 


It is important to note that this widely separated microscopic, mesoscopic, and macroscopic 
spatial scales are not uniquely defined, but rather defined within an order of magnitude. 


The point P could belong to the interior of one of the two phases or could be located on the 


boundary A, , between them. 


In the example shown in Figure 10 the point belongs to phase a . 


‘4 Tf P and P, are any two points of the biphasic medium, a geometric translation that carries P, to P, 


maps B (P) into B (P, ) . In this case it is said that B(P, ) and B(P) are congruent by translation without 


rotation. 
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% 
Figure 10. Sketch of a cubic representative volume element B(P) corresponding to a point P in phase @ 


with position vector X . The position vector of any point Q is ¥+ & 


Figure 10 sketches a cubic REV corresponding to a point P. The position vector of the centroid 
of this REV is represented by X . As consequence, B (P ) can be also written B(x ) . 


The position vector of the origin of the orthogonal cartesian coordinate system shown in the figure 


is 0, and B (0) is the REV whose centroid is this origin. 
Applying the translation given by x, B (0) is transformed into B (x) : 


§ B(0)H3(E)=F + Fe B(z) (3) 


The REV can be expressed as the union of three sets, two of them corresponding to the phases 


and a third one A, B (P ) to the boundary between these phases: 
B(P)=B,(P)VA,,(P)UB,(P) (4) 


We define B i (P ) (=a, f ) such that its intersection with A, B (Pe ) is empty. 


An indicator function can be defined for each phase: 


_| lif Q<B,(P) 
LO\o4 pene (5) 


The volume fraction of phase ! in B (P ) is: 


TET anol =D (P) (6) 
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The volume V (B (P )) = V is independent of the considered REV. However, if the medium is 
non-homogeneous, V (B, ( P)) =V,, ( P) will vary from point to point. 


Because the volume of A, B (P ) is zero, V=V, +V, and we have at every point of the biphasic 


medium: 

f.(P)+ fy(P)=1 (7) 
If the medium moves relative to the laboratory frame of reference, the phase indicator functions 
will be functions of time. 
Let us suppose that the movement of the medium develops in a time scale of a numerical order of 
magnitude greater than the numerical order of magnitude of the time scale of the macroscopic 
description of transport processes in the medium. 
Taking this last time scale as the scale of reference, it is possible to apply the first principle to link 
time scales of different orders of magnitude: the variables belonging to processes with time scales 
at least an order of magnitude greater than the reference time scale, can be considered as frozen. 
Then, it is possible to study the dynamics in the medium using time-independent phase indicator 
functions. This is assumed from now on. 


Let y(t,Q)=w/(t, ¥) be a scalar field in the medium. This is called a local point field. 


From now on we suppose that the local point fields are smooth enough in the interior of the 
phases where they are defined." 


The extrinsic phase average of the field in the point P is defined relative to the whole REV 
volume, being y, (1,0) = y(t,Q) ti, (Q) : 


1 
v.N(P)=4f,,¥(60)1,(0) av 5 
If w (t, P ) =l,y, (t, P ) =f, (P ) and the extrinsic phase average is the volume fraction Fi : 
1 V 
(1,)(.P)=2f,fu(0)aV =Y= 5, (P) ° 
The intrinsic phase average of the field in the point P is defined this way: 
f 1 
Za) ()=7 Jaan ¥ (2) 1, (Q)av (10) 


Phase averages of a given field, either extrinsic or intrinsic, will be called local averages of the 
fields. 

The relationship between extrinsic and intrinsic phase averages of a local point field at the same 
point of the medium is: 


(v,)(P)=f.(P) (vn) (P) (1) 


From (4) and with the symbols introduced in Figure 10 captions, it follows: 


'5 In general, enough smoothness means that the field is continuous and bounded, with continuous and 
bounded partial derivatives in B,(P)« M=a, B ). If this set had been defined closed, the bounding 


hypotheses would not be necessary. 
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(v,,)(t3)=F [je (3) L,(8) dV (3)= 


1 (12) 
7 Iya (+E) L(+) av (8) 
For a cubic REV the extrinsic volume phase average of the field y,, (t X,Y,Z ) in this cube 


centred in the point (x y, z) is given by the following equation: 


v,)(ia r=)? t [f° v,(tx+E,y+z+o)dé Jana u=a,B (13) 


-£) J_f 
2\° 2 


Figure 11 suggests what could happen to the extrinsic phase average (Vi es y, z) = (Va )(X) 


when the lengths of the cubic REV’s sides vary in a thought experiment. 


4K <4 


/ i. ——. 
/ 
{ 


0 ra a 
Figure 11 An example of a possible variation of a local averaged field (y,,)(%) when the side of the 


averaging cube increases. The effects of microscopic heterogeneity appear to the left of £ »» While the 


effects of macroscopic heterogeneity appear to the right of @, as a slow variation of the average when 


increases. 


If € is too small and the point P does not belong to phase @ the average is zero. 


Then, as ( increases the phase average first increases (or decreases) from zero, and then oscillates 


in the microscopic scale of length d , until an intermediate interval of lengths, shown in Figure 11 


centred in € x» 1S attained. 
When £ is not too great or not too small compared to ¢,, the phase average (Vv, )( x) remains 


almost constant when @ increases or decreases. 


If £ increases even more, (v.)(X) begins to vary relatively slowly, now in a length scale L. 


When d « (, <L wecan choose one of these lengths near enough to C, asa mesoscopic length 


scale and L as a macroscopic one. 
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If we average the average (v,)(t x ), we obtain the original average plus one term of the same 


2 
numerical order of (+) v(t.) [60]. 


* 


This last term can be neglected if restriction (R3) is verified, that is if <1 


In this case the average of the average can be put equal to the average: 


((v)) =v) (14) 


The original microscopic field Y can be decomposed into an averaged term (v ) varying slowly 


in the macroscopic length scale L, and a term y'(not necessarily small relative to (v)) 


fluctuating in the microscopic scale of length d and having zero average: 
y=(y)+y' (15 a) 
(v')=(v)-((w)) =0 (15) 


In the literature related to volume averaging applied to transport processes in porous media, the 
decomposition (15 a) is often called a Gray’s decomposition. 


Analogous concepts can be introduced for time averages of fields. 


A time average of the signal g (1) can be defined by z(t)= call, { g(t')dt’ 
t—At 


The time average of this time average is z(t)= 7 ; { | i g(t") irae 
t t-At 


If f) is the time scale of g (1 ) , then z (1 ) is the original average g (r) plus one term of the same 


t'-At' 


0 0 


numerical order of “Jew: Xe) g(t)= a(t) when <I [51]. 
t 


Since some fields are determined referring to surfaces (in general all flows and certain scalar 
fields such as the transmembrane potential) it is necessary to relate the averages with respect to 
volumes with the averages with respect to areas. 


Figure 12 shows a sketch of a REV centred in a point with orthogonal cartesian coordinates 


(x,y,z). 
The REV can be cut transversely with respect to each of the coordinate axes, thus originating 
three families of cutting planes. Those planes belonging to the same family are parallel to each 
other. 
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My 
Figure 12. Left: A plate of thickness dx cut orthogonal to the x-axis in a REV centred on the point 


(x, y, z) . Right: Orthogonal cut to the x axis located atx+& The points of this surface are characterized 


by their coordinates 77 and ¢ . 


The area averaged field in a cross section orthogonal to the axis Z is, by definition: 
be abe 
Vue (Bie) = sfijiv. (t,.xt+é,y+n,z)dndé (16) 
The volume average is related i ne area average through the following equation: 
'p 
(v,)(t y.2)= [7 Paes (om y.ztc)de (17) 
«2 


The line average of the area average is equal to the volume average (known in quantitative 
morphology by Achille Delesse’s rule) [29] [57] [59]. 


For a smooth field y , (t a3 ) the volume average and the area average are related by [52] [60]: 


eee : = ve 
7 LiPo (4x+6.9,2)d6= (tx, 2) +O (t, x,Y,2) <] (18) 


/ 2, 
Here 0 a Ce x,y, z) a denotes the order of magnitude. Analogous expressions are 


obtained by relating the volume average to the other two area averages. 


2 
Then, unless a relative error 0} fs | : 
L 
(W(t. 902) =yacey (t5 25902) = Wyre (6% WZ) =Vry ys (t.25 922) (19) 
The same approximate equality is obtained for the volume and area averages of a smooth vector 
I(t, x,y, Z) or tensor T (t Xx, Y,Z ) field, applying (19) to each scalar cartesian component of the 
vector or tensor field. 
When (19) is applied to the indicator function / 3 (x) of a phase, we obtain the approximate 


equality between the volume and the area fractions of a given phase: 
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Ful P)= fuye(P)= fuse P) = fuey(P) (20) 


Now, let us return to the problem of defining the mesoscopic spatial scale. In Figure 11 a thought 
experiment is suggested for the relationship between the average of a microscopic field on a REV 
and a characteristic dimension of that REV. In general, this thought experiment cannot give rise 
to areal experiment in which microscopic fields can be measured well enough to perform volume 
averages. 


1 if yeB, (x) 
0 if y¢B,(x) 
Contemporary techniques of quantitative morphology based in modern numerical algorithms and 
computer hardware can be applied to determine this function with enough detail to perform the 


However, there is an exception: the indicator function of a phase: [ , ( y) = | 


required volume averages. The details of this nontrivial procedure will not be discussed in the 
present work. 

The arguments raised in relation to the unspecified averaged field shown in Figure 11 can be 
applied, with the intention of determining the mesoscopic scale, to the relationship between 


( i )(2) and the characteristic length of the REV centered at the X point. 


Let us consider now a few useful theorems related with volume averaging of local point fields. 
For this purpose it is necessary to extend the definition of these fields from points interior to the 


phase to points located in the interface A, B (a) . The interface is assumed to be a smooth 
enough surface: fractals are not considered here. 
For any point 5 € A, B (x) , the extension (from phase a ) of the local point field to this point in 
the phase boundary is defined by the following limit: 

v..(t,5) =limy (7, ¥) 

yea 

An analogous definition applies to a scalar field in phase f and to the cartesian components of a 
vector or a tensor field in both phases. 
If Ay, (5) represent a unit vector field defined in points of A, B (x)= Ape (X) and directed from 
phase @ to phase £ we have these formulas relating the average of the gradient with the gradient 


of the average, and the average of the divergence with the divergence of the average [15] [60]: 


(Vy, )(t.%)=V(w,)(t.%)+ W.,(t,5) fi,,(S8)dA(S) (21) 


1 
TB) lento 
(7 F)6)= VME) TT In 


Interchanging a@ and f, formulas for (Vv V,)(t,¥) and (V-Fg)(t.¥) follow from (21) and 


1 


n 


ap (5)-Ja(t5)dA(S) (22) 


(22), with #,,(5)+i,,(5)=0 VS eA,, (2) 
Representing by § ( Aes (x)) the area of the interface between the phases in the REV B (x) , the 
surface integrals that appear in the right side of equations (21) and (22) can be rewritten as a 
product of an area fraction Hos ( x) = 5(Aw()) by an area average 


V(B(x)) 
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Then, the equation for the divergence of a vector field (22) in this equivalent form is: 

(V-7,)(t.¥)=V- (in (t.%)+ Xap (¥) J, (t,¥) (23 a) 
In (23 a) appears the area average of the current density that crosses the interface from phase a 
to phase #. As consequence we have: 


Ha it Fiap (8) J (tS) 4 A(5) (24) 


(V-FN(42)=V-(p)(6¥)- Za 0(8)4.(63) 036) 


If w, gtepresents the velocity field of the interface relative to the laboratory frame of reference: 


(Se Net)= 2 We) O8) PJ gq MelO8) ap 8) Raol8) 4A(8) 


In the time scales considered here it is possible to put w,, = 0 everywhere. So, the volume 


I AEx)= 


average of the partial derivative respect of time is equal to the partial derivative of the average: 


(Su. \(e8)-ZaNea) (5) 


To formulate the next equation, we begin by a decomposition (equation (15 a)) of a scalar field 
and a vector field into an averaged term varying slowly in the macroscopic length scale, and a 
term (not necessarily small relative to the averaged field) with small scale fluctuations and having 
zero average: 


(v= Wa) *¥% Hiei) t, 
Then we have: 

Vs, i=.) (in) + Vi iu) (26) 
Finally, we have the following lemma, being M an operator that maps a vector space into itself: 


Lemma: If for V b the operator M verifies b-M (a ) =M (6 ) -G , then M isalinear operator. 


These results can be applied to the volume average of the microscopic fields and equations that 
describe the transport processes of a substance-like quantity (mass, charge, thermal energy, etc.) 


of volumetric density g ( y) : 


If the current density of the quantity is represented by i(t ; y) , the balance law is: 


x HO . 
V-F(t9)+ 2 a(69)=0 (27) 


Now we introduce a Fick-like law in terms of a generalized conductivity tensor K (¥). 


symmetric and definite positive, and a generalized scalar potential w(t, ¥) (concentration, 
electric potential, temperature, etc.): 

i(t,¥)=—-K(¥)-Vy(t,5) (28 a) 
The generalized formulation of classical transport laws given in equation (28 a), assume an 
instantaneous response of the current density to the applied scalar potential gradient. However, it 
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takes some time to the currents to reach the values predicted by the classical transport laws. This 
can be important when high enough frequencies and short enough wavelengths are involved. 
A possible modification of (28) is to introduce a relaxation time 7 of the current to the quasi- 
steady state given by the classical transport law [25]. Then, the new constitutive equation is: 


rE 5(1,5)+ I (0,3)=-R (3) Vo (63) (28 b) 


If the medium is isotropic, i represents the unit tensor, and K ( y) is a local scalar conductivity: 

K(y)=K(y¥)I (29) 
Possible cross effects due to the transport of other quantities are, for the moment, neglected. 
However, they will be considered in section (E). 


) 
qd (t, y) andy (t, y) are related by a generalized volumetric capacity field C = AG q: 
y 


O O 
Cae ere eee 
ay At y) OAV y) (30) 


From (27), (28) and (30) the following diffusion-like (that is, parabolic) equation is obtained: 
=\ 0 % sf e 
C3) Sw(65)=V-(K(5)-Vy (4.5) (1a) 


If the relaxation of the current is considered, from (28 b) and (30) we obtain this hyperbolic 
equation (which is one of the forms of the so-called telegraphy equation) [61]: 


CO Svs) Sv (3) ]=v(RO)VH3)) (31 b) 


These considerations have a certain interest in our case, because recently a modification of the 
bidomain models of electro cardiology, was proposed modifying Ohm’s law by the introduction 
of a relaxation time for the electric current density [45]. 


Nevertheless, from now on we work with the classical version (28 a) of the generalized transport 
law. So, let us apply equation (31 a) to phase of a biphasic medium. 


Introducing a representative value C of the generalized capacity: C ( y) =C “0 C ( y) 
If K »o 1S a representative value of the generalized conductivity: K ( y) =K “0 K ( y ) 


Both C ( y ) and K ( y ) are dimensionless. Considering the spatial scale d , and changing the 


spatial variables to dimensionless ones by y =d i y, so that V =—V, (31 a) can be rewritten: 
il 
syd g = - 
tame €(3)—w(t.3)=¥-(&(5)-Yy (1.5) Gc) 
dC 
ty, mic — a (32) 
, K 
7a) 
Here?,, ,,-18 a characteristic time of the equalization process described by equations (31), and 


0: 8 . ‘ . ar 
— is a representative equalization coefficient. 


yay) 
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Suppose that ¢ is much less than the macroscopic time scales of the averaged fields '°, and 


b,mic 
the description is done relative to this last time scale. Then, the second principle to link time scales 
can be applied. The variables belonging to the process that evolves with a time scale at least an 
order of magnitude smaller than the macroscopic variables of reference, after a short transient, 
can be considered as relaxed to equilibrium with these variables: the microscopic process can be 
described as in quasi-steady state. 

So (25), written for phase /, is simplified to: 


V- jy (t,¥)=0 (33) 
Applying (33) to phases a and it follows, considering (23 a): 
(V-..)(t.%)=0 (34 a) 
(V-Jp)(t-¥)=0 (34 b) 
V-(Fa) (0%) = Xap (2) J. (6%) (35 a) 
V-(is)(t4) = Xap (%) J, (16%) (35 b) 


As will be seen in more detail in section (F), when charge transport is assumed to be quasi-steady 


at the microscopic scale, J, is the microscopic electric current density in the intracellular space, 
cP in the extracellular space, (i) | , and (Fay =J p are their extrinsic volume averages over 
the REV B (xX ) of centre X, y, B (x) is the local fraction 7,,, (x) of excitable membranes, and 


w. (Aa ) is the macroscopic density of membrane current MA m> then (35 a) (35 b) are identical 


with the equations (B) (5) of the bidomain models of electro cardiology. 


The last problem to be considered in the present section is the so-called closure problem: to get 
an approximate solution, based on the separation of scales, of the coupled micro and macro-scale 
transport problem for a given phase. The goal is to obtain the following macroscopic averaged 
form of the microscopic, generalized Fick’s law, in terms of macroscopic variables only: 


= wy ~ Ee f Sy 
GH Nt2)=—F, Kv (%)-V(y,,) (t,x) (36) 
In (36), K neff (x ) is an effective macroscopic property resulting from the microscopic property 


K ii (x ) after applying the closure operation [7] [26]. Even if the phase behaves as isotropic at 
the microscopic scale, when the averaging process is applied, the resultant effective property 


K weg May be a tensor quantity if there is a degree of anisotropy related with phase geometry. 


Because the macroscopic flow ( ie ) is an extrinsic phase average, while the macroscopic force 


f 
V (y ) is an intrinsic phase average, and because from approximate equation (20) it is possible 


to work with volume averages instead of area averages when dealing with current densities and 


av.) (63) 


being Li, (y) 4 time scale of variation of the 
Ot 


16 : 
For example, if ee | (v,)! ial/ 


i 
averaged field (y ‘a (t JX ) , defined at a given time and REV location. 
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the phase fractions associated with these currents, then it is reasonable to have the phase volume 


fraction f , as multiplier of the effective mesoscopic conductivity. 


There are two main approaches to the deduction of (36) for the transport of a substance-like 
quantity. 

A first, less powerful one, is based on a reciprocity theorem and a very general closure hypothesis. 
Two different transport processes of the considered substance like quantity, within phase a . The 


microscopic fields associated with these processes are Y,,, ie and Y,,, a The microscopic 
Fick-like transport law is (equation (28)) vA =-K a VW, Because the tensor K ~ 1S Symmetric: 
VV , Jos = VV ‘K, VW =—VWas ‘K,, Vy =VYn ‘ ce (37) 
Because the microscopic processes can be assumed to be relaxed to equilibrium with respect to 
the slowly time varying macroscopic processes, V - he =0 so V- (Ca ia) ZVW.: Vs 

From this last equation and from (37) we have: 


V-(Yur Jaa) ®V (Wor i) (38) 


From theorem (22) applied to the vector field W,, ie : 
ve Wai ee aye War ie ee War fiz Jad A (39) 
V “Aap 
In (39) the normal component of i at the interface is equal to the current density j,, through 


this interface 7, Blea J a 


We make this assumption A from now on: The interface presents a very high resistance to flow 
relative to the resistance in the interior of phase @ , so ina first approximation any surface integral 


involving 71, B* 7m and bounded multiplier fields like Y, can be neglected. 
Then (V ; (Vin ja) aV- (Vo in) and (Vv (War ju) aV- (Vr ia) » SO from (38): 


Ve (Ver is) aV- (Var ia) (40) 


From theorem (26), (Var ia) i (va) ce ) au (v.1 ee 


Because V - (ia =O (assumption A): V- (Wea i) —V (yw...) ‘(Fu2) +V.- (va! oe 


d 
The term V - CA >) can be neglected because it is proportional to on <1 []. 


* 


So V-(War Jar) V (War) (Gar) and V(r Jn) © V oa) * Gar) 


Then, from these last relationships and from (39) the following reciprocity theorem is obtained: 


fo faz fe oft 
V (War) Gar) ®V (War) * Gas) (41) 
Now, we make a general closure hypothesis: the extrinsic average of the current density field is 
an unknown function of the gradient of the intrinsic average of the potential field. 


j.=M[Vv.)" | (42) 
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From (41) and (42), considering that the vectors V (Wor i and V (yo)! can be any vector in 


. f 4% f f 4% f 
physical space: Vv (War) -M ky (War) ] xV (War) -M ki (Wat) ] 
Finally, applying the Lemma it follows that M must be a linear operator: 
Ja =—ta Keg Vay (43) 


In (43), K aef 18 an effective macroscopic property resulting from the microscopic property K, 
Averaging (37): (Vv ° Ten) + -(Vv War ‘ K, VW) a -(v Wao ; K, VW) = (Vv Wo : ia) 
In particular: (Vy, ja) =-(V Vv, -K, ‘VyW,) 

; ot if ope 
From these last two equations, from V- (Vor jus) aV a) ; (i,2) and the analogous 


cs fv s 
equation V- (Ver in) x V(W.) : ( dae from the positive definite nature and symmetry of 


the microscopic conductivity tensor, it follows that the macroscopic conductivity tensor is also 
definite positive and symmetric. The deduction of the Fick-like macroscopic law is complete. 

If scale and regularity constraints on geometry and fields are applicable, the result is applicable 
for any microstructure of the continuum. The components of the effective conductivity tensor 


K., must be estimated from the results of experiments carried out on a macroscopic scale. 
However, the above demonstration is not constructive: we don’t have a procedure to construct 


K, a {rom the microscopic tensor K,, and the geometry of the regions corresponding to phase a 


. As consequence, the components of the effective conductivity tensor K,, a must be estimated 


from the results of experiments carried out on a macroscopic scale. This limitation is removed by 
the second closure method. 


This second approach introduces a representative unit element (RUM) centred at X , assumes that 
Bs ne ae es 
y (1 ; y) = b,, ( y) -V (y a (1 1X ) there, and deduces a non-homogeneous elliptic boundary 


value problem for the unknown field b, (y) to obtain the desired upscale leading to the 


macroscopic version of generalized Fick's law. 
The RUM is a volume element with, in most cases, spatial scale @, (it can be a cube with side 


length equal to the mesoscopic scale) but its contents don’t necessarily coincide with the contents 
of a REV (although they may coincide). 

A RUM is constructed from images and conductivity measurements of the biphasic medium and 
is intended to capture the most important geometric and physical properties of the complex 
medium. When this medium is composed of regions disordered at the microscopic scale but 
relatively homogeneous at the macroscopic scale, for each homogeneous region a RUM may be 
composed from an averaging of the region random properties, geometrical and physical. If the 


statistical properties of the biphasic medium change, in principle a RUM family { B, (%)} 
(centred at different locations X ) can be constructed. 
From now on spatial and area phase averages will be taken over any RUM B, (x i. Because 


restrictions (R1), (R2) and (R3) remain valid for a RUM, all the theorems and bounds obtained 
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previously for REMs can be applied to RUMs. In particular, the microscale potentials and flows 
can be taken as if they were relaxed to equilibrium with the macroscopic, averaged, variables. 


From decomposition (15 a) applied to the phasic potentialy, (e y)in the RUM, from the 


microscopic Fick-like law (28), and from the quasi-steady equation (33) for the microscopic flow 
density it follows: 


V-(K. (5): Vy. (4,5))=-¥-(K, (8)-V v2)’ (63)) (44) 
This microscopic balance equation shows a coupling between the microscale field y/ (t, y ) and 
the macroscale field (y,, \ (t,x)for every Ye B,, (%). Here By, (X) is the region occupied 


by phase @ inside the RUM B, (x) 
We make assumption A made above (the interface presents a very high resistance to flow relative 
to the resistance in the interior of phase @ , so ina first approximation 7, B° hi =0). 
Then, we have this boundary conditions at the interface inside the RUM: 

fi, K,,(3)-V iy (t,3) =i, Ky (3) (Vw) (6%) VIE Ag (X) (45) 
This microscopic boundary condition shows a coupling between the microscale field y/, ie y) 


and the macroscale field (y,, \ (t,x), now forevery ye A,, (x). 


In the non-homogeneous microscale elliptic equation (44) and in the non-homogeneous 


microscale boundary condition (45), two source terms appear, both involving the macroscopic 
averaged field (y, \! 
Now, we add this boundary condition in the points belonging simultaneously to phase @ and to 


the surface 0B, (x ) that separates the RUM from its environment: 


Va (t,9)=F,(t,9) VIE OBy (%)ABy, (¥) (46) 
For developing an effective closure, we assume that (46) can be substituted by periodic boundary 
conditions along three mutually orthogonal directions. 
The nature of the balance equation and the boundary conditions at the interface suggest the 
introduction of the following tentative solution (ansatz): 


y,| (t,3)=b, (5)-Viw,) (8)  — V¥ E Ba, (X) (47) 
The next step is to introduce the decomposition K, (¥) = (K, . (x) + Ke (¥) as well as the 
ansatz (47) in (44) and (45). 
After some formal calculations, considering that (K * i and (Wa ‘behave as if they were 


constant in the RUM and its neighbourhood, we obtain the following non-homogeneous elliptic 


equation and boundary conditions in terms of the unknown tensor field Vb, ( y) 5 


V-(K, (5): Vb, (5))=-¥-(K, (5)] V¥ € By, (X) (48) 

fi,p°K,(3)-Vb, (¥) =i, 9K, (5) VyEA,,(X) (49) 
Instead of the boundary condition (46) we introduce the periodic boundary conditions: 

b,(F+h,)=b,(F) V5 eOB, (%) ABs (¥) (50) 
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From (48), (49) and (50) the field b, (¥) can be determined unless a constant vector. 
Taking the average in both members of the ansatz (v7. (t, 5) = (b, (¥)) .V (vy, ‘a (t, x) 
However, we know that the average (v.! (t, 5)) =0 and VW, 7 ( 1.) may be equal to any 
vector, so we must add a further restriction to the field b, ( y ) : 

(b, )=0 (51) 
The non-homogeneous elliptic boundary value problem (48)-(51) can be solved by means of 


computing packages to obtain a numerical solution for b, (¥ ) in the RUM. 


To deduce the macroscopic version of the generalized Fick’s like law, we return to its microscopic 


version I, =-K, ‘Vy, Because Vy, =(Vy, \/ +Vy, and Vy,’ =Vb, (Vy, Sy we 


have: 
iu (1,3) =-K, (¥)-(F+¥b,(¥))- (Wye) (1.3) (52) 
Taking averages over the RUM in both members of (52) we obtain: 
(Fa )(t-¥) = fa (%) Keer (2) (V vey’ (44) (3) 
The effective macroscale conductivity K a eff (x) verifies: 
fa (®)Roseg (8) =(K. (3)-(F+ Vb, (5)))(%) (54) 


The closure operation is now complete, the macroscale Fick’s like law is substantiated and the 
effective generalized conductivity tensor is expressed in terms of microscale fields, either 


assumed known, as is the case of K ns ( y ), or determined solving an elliptic partial differential 


equation with suitable boundary conditions, as is the case of the vector field b, (¥). Analogous 


formulas can be derived for phase /. 


To apply the results obtained in this section, we must discuss first the electromagnetic foundations 
of the bidomain models of electro cardiology. 


(E) Microscopic’ electric flow in cardiac tissue and microscopic* boundary conditions 


A quasi-stationary version of Maxwell equations, suitable to describe the microscopic (in the 


sense of the d -scale) processes in cardiac tissue, can be obtained introducing simplifications to 
the equations of the classical electromagnetic theory of continuous media. 

These equations are obtained applying an averaging process to the Maxwell equations formulated 
at the atomic scale. 

To consider this additional averaging process we must introduce a fourth scale, atomic or a - 
scale, and a fifth micro averaging scale or a” - scale, in addition to the three scales already 


introduced previously: the microscale or d - scale (microscopic in the sense of Histology), the 
mesoscale or ¢,- scale, and the macroscale or L - scale (reaching the lower limit of the 


macroscopic size in the sense of Anatomy). 
In Physics, the version of Maxwell's equations, in thea-scale, is called microscopic. To 
distinguish it from microscopic in the histological sense, we will call it microscopic’. 
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The resulting averaged equations in the a” - scale, instead of macroscopic as they are called in 
« 17 


Physics, will be called macroscopic’. 
To smooth possible abrupt variations when whole molecules exit or enters the small averaging 
volume, the averaging is done taking the convolution of a local property 77 (t aa ) , defined at 
the a -scale in the region occupied by the biological material, with an isotropic normalized 


weighting function w(F ) : 
(7), 03)= [J o—sDnler Fav") 


The weighting function tends gradually to zero as the distance from its centroid increases. 
A possible example is: 


my 2 
Ze, ts 1 _ 
w((\7 5\I) = exp |7 || 


The weighting function must be smooth enough to enable a development in a rapidly converging 
Taylor series over distances in the a” - scale where the weighting function differs significantly 
from zero (in the example of weighting function the linear dimensions R ~ a” is of order of 
10 nm in condensed matter like biological tissues) [20] [48]. 


As a thumb rule, if ais a characteristic length of the atomic scale and d is our microscale, then 


in a first approximation a’  ./ad 


If azO.lnm and d~1000nm=1 wm then a ~l10nm. A microscopic scale (in the 
histology sense) of linear dimension d ~ 1 zm could be considered a reasonable average estimate 
of cytosol distances between neighbouring membranes or between membrane less organelles in 
the cytoplasm, or between collagen fibres and cell membranes bordering the extracellular space, 


given the geometric details of the intracellular and the extracellular spaces in cardiac tissue [1] 
[44]. 


To take the convolution the field 77 (t, ar ) must be extended beyond the region occupied by the 


biological material (a region blurred in the level of resolution corresponding to the a - scale) to 


the whole R ° , but this is not a big issue due to the fast-decaying behaviour of the weighting 
function [48]. 


The Maxwell microscopic* equations are given in terms of the electric field é (t, r ) , the magnetic 


induction field b G r - the total local charge n(t, r) and the total local electric current density 
I(t r ) . The local source free equations are (V A (...) represents the curl and V “(s) the 
divergence): 

VAg(n?)=-Sb(17) V-b(t,F)=0 


The local source related equations are: 


'7 The method of averaging considered here was developed mainly by Roussakoff in his Ph D Thesis and 


is based on a modified version of spatial averaging. The a” -scale averaging of the microscopic* Maxwell 
equations begun with Lorentz work at the end of the 19 century and continued, during the 20 and 21 
centuries, with increasing rigour. It incorporates quantum mechanical systems at the atomic scale and does 
not necessarily imply an explicit time or ensemble averaging [48]. A complementary approach was 
developed by Robinson [43]. The 1998 edition of the book by Jackson summarizes both approaches [20]. 
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ot 
In a classical approach to microscopic* electromagnetism the charges are point like, so the fields 
n(t.7)= > 4,6(F 7, (t)) and j(t,7)=) q,¥, (1)5(F -F, (t)) are given in terms of Dirac’s 
j j 


&V-8(t,7)=n(t,7) Vad(t,F) =n (oZetur)i7)] 


distributions. 
Averaging the microscopic* equations, interchanging the order of averaging and taking 


derivatives of the fields, and introducing E(t, y) = (é), (t, y) and Bit, y) = (6) (t, y) the 


following macroscopic” (d -scale) version follows: 


VE (t,3)=-< B13) (1) 
V-B(t,¥)=0 (2) 
&V-E(t,¥)=(n), (1.5) 3) 
VA B(t, ¥) = My [« E(r, +i) (t, 2) (4) 


These equations are defined in the region (2 occupied by the biological tissues. This region is 
assumed to be sharply defined in the d -scale. 
The internal sources of the fields are related with the average charge density (7), (t ; y) and the 


average current density ( J ) (t : y) and their external sources can be related with stimulating 
a 


electrodes (perhaps located at the boundary of the region Q), stimulating coils (located external 


to Q and perhaps in its neighbourhood), NMR equipment, nearby power transmission lines or 
electromagnetic fields from antennas. 


Due to the selection of the a’ -scale, the resultant macroscopic’ averages (7), (t, ¥), 


> 


a (1, y) ' E(t, y)= (é) (t,¥) and B(t, y) - (6), (t, y) vary significantly only in the d - 


scale, not in the a‘ -scale. 


In biological tissues almost all charges are linked in atomic structures. Free electrons, as can be 
found in crystals, and free nuclei as can be found in plasmas, are in general not found. '® 

When the weighting function is developed in a Taylor series, the averaged charge and current 
fields are expressed as sum of multipolar terms: 


(7), (1,3) = p(t. 3)-V-P(t,9)+V-O(t, 9) +... 


+ 4 S50 Op. ss fp 
(i). (t,¥) =] (1,5) +— PUL y)+V AM (t, 9) +... 
If the atomic structures (atoms, molecules, macromolecular structures) have a nonzero net charge, 
then they contribute to the macroscopic’ free charge density Plt y) as particles. They contribute 


to the macroscopic polarization P (t, ¥) and to the polarization (bounded) charge density 


'8 Proton jumping and redox transport of circumstantially free electrons are possible, but unimportant 
exceptions for our purposes. 
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V-P (t, y) as structures composed of several particles. The tensor O (t, y) of electric quadrupole 


moments and higher order terms usually can be neglected [20]. 
Defining the electric displacement field, the averaged equation (3) can be written as usual: 


D(t,¥) =e E(t, ¥)+ P(t, ¥) (5) 


V-D(t.3)=p(t.3) © 


In the multipolar development of 7 i) (t, y) the term I(t, y) gives the current density of free 


moving charges or conduction current density. The next term ar P (¢ , y) gives the polarization 


current density due to internal movements of charges inside each atomic structure and to the 
movements of the centres of charge of this structures when they are polarized. The term 


VAM (t, y ) is the magnetization current expressed by the curl of the magnetization field 
M (t, y) . Higher order terms in the development of (7) (t, ») are usually neglected [20]. 
Defining the magnetic field, the averaged equation (4) can be written as usual: 


os I eye, tis 
Bley es ee) ee) (7) 


0 
a ee ee ee ee ee 
VAH (4.5) =i (t.9)+— DUty) (8) 
Equations (1), (2), (5), (6), (7) and (8) are the version of Maxwell equations due to Hermann 
Minkowski. 
The current density I(t a) ) is the flux vector of the density p(t : y) of free charges. !° 
From (6) and (8) we derive: 


SA ue 0 = 
V-F(45)+ 2 (4 y)=0 (9) 


The scalar field pit, y) and the three vector fields P(t,y), i(t,¥) and M (t, y) specify the 


electromagnetic properties of the biological tissue and are related with its mechanical and 
chemical properties. 


From an electromagnetic point of view, p(t.y), I(t, y), P(t, y), M(t, y), E(t, y) and 
B(t, y) (and, of course, the derived fields H(t,¥)a nd D D(t, 


the electromagnetic state of the biological tissue. 


y) ) are functional fields that give 


Let us summarize the Maxwell equations in the version due to Minkowski, before continuing with 
the formulation of the equations of the electrodynamics of continuous media for a biological 


tissue. In terms of the fields E(t, y) ,B (t, y) ,D (t, y) and H (t, y) we have two homogenous 


equations (equations (1) and (2)) and two non-homogeneous equations (equations (6) and (8)): 


ae a ete . 
VAE(t.¥)=—= B(LS) V-B(t, ¥) =0 


'9 Tn the sense that the rate at which the quantity of free charge that crosses per unit time an element of area 
dAis 7-dA 
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- ee eee ee ee 
V-D(t,¥)= p(t. ¥) VAH(t9)= J (49)+— DUty) 


or to 0 a 
Besides, we have the equation (9) of conservation of free charge: V: J (¢ ; y) ze at Pp (1 +y ) =0 


From these equations the following boundary conditions at a smooth interface OQ.,, between 
two regions Q), and O; are derived [16] [23] [47]. The sub index ¢ represents the tangential 


components of the fields relative to the interface, and n represents its normal components, while 


@, is the density of free charges at the interface and I, (t , y) is the density of pellicular currents 


there. We have, V y€0Q,,: 


Ej, (by) =En(hy) (10a) 

B,,(t,9)=B,.(t,5) (10b) 

D,, (t, ¥)—D,, (t.¥) = @, (t, ¥) (10c) 

H,,(t,9)-H,.(t,9)=3,(t,9) (10d) 
eM RRS ie ee coe : 

Jn (t.9)— jo (t,9) = 5% (1, ¥) (10e) 


To these general boundary conditions are added equations (5) and (7) between the derived fields 


and the main fields, which involve the polarization P (e y) and the magnetization M (t, y ) fields. 


From now on we will neglect the magnetization field because it is not necessary to the formulation 
of bidomain models in electro cardiology considered in the present research. 


To continue, we must relate the fields P(t, y ) and I(t cy ) with the properties of the tissue and 


the electromagnetic field there. 
For the conduction current density, we assume the following decomposition: 


i(t,9)=6| E |(t, 9) +7.(63) (11) 
The first term on the right-hand side of the equality is a generalized version of Ohm’s law, written 
for a linear, dispersive, anisotropic, and non-homogeneous medium. 
The causal linear tensor operator 6 [...] is a generalized electric conductivity for a linear, 
dispersive, anisotropic, non-homogeneous medium, that approximates part of the conductive 
properties of the biological tissue []. The components of the vector field 6|E ] (e y ) are: 


A — 3 y 
6| E| (t, y) = by | C7; (u, ¥)E, (t—u, y)du i=1,2,3 (12) 


J=l co 
The second term ie (t , y) encompasses the bulk effect of the concentration gradients of mobile 


ions in the fluid spaces delimited by the cytoskeleton and the micro-trabecular network, the 
entrainment due to the possible bulk flow of tissue fluid and localized ionic flows through open 
channels in cell membranes, driven by electrochemical potential differences and possibly related 
with nonequilibrium phase transitions in biological materials [14] [17] [41]. 

For the polarization field we assume the following decomposition: 


P(t,3)=6,4[E|(t.3)+ B(t3) (13) 
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The causal linear tensor operator x [. ; | is a local generalized dielectric susceptibility for a linear, 


dispersive, anisotropic, non-homogeneous medium, that approximates part of the dielectric 
properties of biological tissues. The components of the vector field y [ E ] (t, y) are: 


A | 3 : 
ZEL(LH=>) | nH (t-uj)du  i=1,2,3 (14) 


j=l 
The additional term P(t, y) encompasses bulk effects and localized membrane effects on 


bounded charges related with conformational changes and other cooperative effects in 
macromolecular structures [12] [14] [40]. 


Introducing the dielectric permittivity operator € [...] =& (i & +7 [..]) equation (5) may be 


written: V-(2[E])(.5)= p(t.5)-V-B (5) 
From (11) and D(t, y) = [ E |(t. y) + P. (t, y) , the Ampere-Maxwell law (equation (8)), can 


be written as follows: 


ey! ANC in OO) ea: elie, Se x ae: se 
VAA(t3)=H(45)+ OP (t5) +e E|(t5)+ é| E |(1,¥) (15) 


Ot 


The operator fields 6[...] é[...] and the magnetic permeability operator field ii [eal (when 


tissue’s magnetization is important) are structural fields: they give the electric and magnetic 
properties of the medium, which in this case is assumed to be linear, stationary, local, dispersive, 
non-homogeneous and anisotropic. 

However, the local electromagnetic description is not complete until the current 


ea eee 
Js (t ; y) + a Ved (t ; y) is specified. In the usual formulation of the approach known as the volume 


conductor theory this term is included in the designation impressed current or primary current, 
depending on the level of resolution. 
When cellular or sub cellular level of details are left without resolution, everything in 


i(t.y¥)= 6| E \(t, ¥)+ j,(t, ¥) not included in G| E(t, y) is considered as primary current 
[17]. 


A further distinction is done between currents through active biological membranes themselves 
and secondary currents related with interfaces [36] [37]. 

When the purpose of the theoretical description is to enable the identification and quantification 
of the spontaneous current sources in the tissues, from electric or magnetic measurements, these 
currents are considered as given by a parametric model whose parameters must be estimated from 
experimental data [16] [17]. 

But when the purpose of the theoretical approach is to describe the genesis of action potentials in 
excitable tissues, due to electric or magnetic stimulation, a mathematical model that relates these 
currents with certain cellular and subcellular variables, must be introduced [2]. 


From a macroscopic point of view, the electrical bulk properties of a tissue are measured as a 
function of the frequency of an alternating current signal (alternating current spectroscopy). 

If the amplitude of the input signal is small enough, the directly measured quantities (electrical 
conductance and capacitance) are independent of the signal amplitude and depend only on the 
applied frequency. 
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To avoid the effects due to the polarization of the electrodes an alternating current of controlled 
amplitude and frequency is generally used as the input signal. 

The amplitude of the current is selected so as not to stimulate excitable cells that could be part of 
the sample. 

However, if the volume fraction occupied by excitable cells is negligible, the non-linearities 
associated with the appearance of action potentials do not influence the result of the measurement. 


A four-electrode method is used to measure the global conductivity of the tissue. 
Two of them inject current (at a distance d from each other) and the other two measure potential 


difference (they are at a distance L from each other and present a negligible drain of current). 


S 
Thus, from the measured conductance G(a) and from the relation G(o) = raACOE the 
conductivity 0 (a) equivalent to a homogeneous portion of tissue is obtained. 
S 
Analogously, from the relationship C (a) = 7 EE, (a) the equivalent relative dielectric 


permittivity €, (a) can be determined. 


Figure 13 left suggests how these measurements are done in a sample of biological tissue placed 
between two flat electrodes. 
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Figure 13. Left: A sample of biological tissue is located between two square electrodes of area S and 
separated by a distance d . Two wire electrodes are shown in the sample, separated by a distance DL. 
Right: Representative measured values of a global relative dielectric permittivity (left ordinates) and a 
global conductivity (right ordinates) as functions of the frequency of oscillation of an applied electric field, 
for a sample of a biological tissue. 


The right side of Figure 13 shows the variation of a global relative dielectric permittivity €, (a) 


and a global conductivity o(@) , in a sample of biological tissue like the one shown in the left 


@ 
of the same Figure 13, as a function of the frequency of oscillation f = ae of an applied external 
a 


field. 


Because we are interested in the electric behaviour of cardiac tissues below 10 kHz, the results 
shown in Figure 13 suggest that, at least in principle, it is possible to neglect dispersive effects in 
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the bulk conductivity. But dispersive effects in this interval of frequencies produce an almost 
three orders of magnitude variation in the dielectric permittivity, so this cannot be neglected. 
Furthermore, we assume that, apart from localized excitable membranes, the rest of the medium 
can be described as electrically linear [12] [14] [16]. 


Then, the current density can be re written introducing a conductivity tensor field G (¥), 


symmetric and definite positive, which remains locally invariant in the time scale of the studied 
processes: 


i(t,¥)=6(3)-E(t,5)+7,(t,¥) (16) 


For non-homogeneous electrolytic conductors, the electric field is not the only cause of the 


conduction current in bulk es ee y ie These additional forces are introduced through another 


vector field E P (t, y ) (called the Planck field) such that [4] [16] [52]: 7° 


foun (#3) = &(3)-(E(t,5)+ E, (1,5) (17) 


From equation (2) it follows that the induction field has a vector potential: 
B(t,3)=VAA(t,5) (18) 


From (1) and (18) it follows the electric field must be given by: 


eyo a Z Ode Ls 
E(t,9)=-Volty)-7 AWS) (19) 


Here g(t : y) is a scalar potential. Under quasi-stationary conditions we regain the electrostatic 
relationship E(t, y) = Volt, y) : 


In any case, when applying formula (17) to the density of conduction current in bulk, in principle 


O55. % 
both terms in the right-hand side of formula (19) must be employed. The term a A(t , y) allows 


us to describe the electric field induced in the tissue by an external time varying magnetic field 
and the term Volt, y) is necessary to verify the boundary conditions in the interface between 


tissue and environment, even if there are no electrodes producing their own electric fields there. 
These matters will be discussed later. 


Now, let us consider the validity of the quasi-stationary (QS) approximation to the Maxwell- 
Minkowski version of the @™ -scale (macroscopic*) averaged electromagnetic equations, applied 
to a biological tissue. 

For a linear, homogeneous, and isotropic volume conductor, quasi-stationary conditions may be 
assumed if propagation effects, capacitive effects and inductive effects produced by bioelectric 
sources distributed in the volume conductor are neglectable, and the normal component of the 
electric field in the boundary of the volume conductor can be set equal to zero [35]. 

These results can be extended to a linear, non-homogeneous, anisotropic, and dispersive volume 
conductor, as will be shown in the Appendix. 


20 A short introduction to the origin of Planck’s fields can be found in the classic book by Abraham M and 
Becker R (1962) The classical theory of electricity and magnetism, Blackie & Son, Glasgow, UK. 


42 


August 2021 


Assuming as a first approximation the validity of the quasi-stationary conditions, the Faraday’s 


induction law given by equation (1) reduces to V AE (t, y) ~0 so there is a scalar potential 
g(t, y) such that: 

E(t,y)=-Vo(t, ¥) (20) 
The Ampere-Maxwell law given by equation (8) reduces to V A H (t, y) ~ I(t, y ) so taking the 


divergence and considering the bulk current density 7! ae (t ; y) given by equation (17) we 


obtain: 
V* Fun (ts 9) = V- (6 (¥): (E(u ¥)+E,(43)))=0 (21a) 
Considering the bulk relation D(t, od é[E E \(t, y) 


(¥-é[2)(.5)=0(.5) a2) 


Some people [49] object to the use of equation (21) when it is written this way, as is often done 
in the volume conductor theory [35] [36]: 


equation (5) gives: 


V-(6(¥)-Voe(t,9))=V-7.(t,5) (21b) 
where Ie (t : y) is assumed to be independently given. For a homogeneous and isotropic 
conducting volume, the scalar potential results to be inversely proportional to the scalar 
conductivity. If this conductivity tends to zero and - (t , y) is independently given, the potential 
should increase without bounds [49]. We see that with (21a) this problem disappears, because 
E P (1, y) remains bounded and so this equation reduces to the trivial equality O=0. 

However, the potential always verifies the equation (v é [Vo] \(e ¥)=—p(t,¥) independently 
of the value of the conductivity. When the medium is a volume conductor, both equations 
(v -é[Vo] \(e ¥)=—p(t,¥) and V-(6(3)-Vol(t, y)) ay (6(5)-E, (t, y)) are compatible and 
simultaneously verified. 

Amongst other feedbacks that appear in a conductive medium, part of the charge density adjusts 
to enable the establishment of a global pattern of current flow in the volume conductor, under the 


influence of the restrictions stemming from boundary conditions and the influence of both, the 
spatial and the temporal distribution internal and external sources of current flow. 


Under certain circumstances, when applying the Maxwell-Minkowski equations to a biological 
tissue, the quasi-stationary approximation must be loosened a bit. 

This is done in two directions: one is called quasi electro stationary (QES) and the other is called 
quasi magnet stationary (QMS). 


In the QES approximation - B(t, y) is neglected in Faraday’s induction law, but the density of 
t 


displacement current - D(t, y) is retained in the Ampere-Maxwell law. 
t 


21 That is, in the cytosol or in the extracellular space. 
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On the contrary, in the QMS approximation = B(t, y) is retained in in Faraday’s induction law 
t 


and - D(t, y) is neglected in the Ampere-Maxwell law. 7 

These two approximations may be more easily justified working in the frequency domain, with 
phasor vector field solutions of Maxwell-Minkowski equations, corresponding to monochromatic 
oscillations [23] [47]. 

The case of a linear, homogeneous, and non-dispersive conductive and dielectric material was 
carefully studied by Haus and Melcher in the time domain [18]. A generalized perturbation 
approach to Maxwell-Minkowski equations appears in [23] and a more detailed version suitable 
to describe biological tissues appears in the Appendix of the present work. 


In the framework of a full QS approximation, with the equations (18), (19) and (20) in terms of 
local point fields, we have the starting point to apply the phase average method and derive the 
bidomain equations of electro cardiology. 

Before doing that, let us analyse more in depth the relation between the bidomain models of 
electro cardiology and the classical electromagnetic theory of continuous media. 

To do this, it is necessary to distinguish two very different levels of averaging. 

One level of averaging allows to pass from the a-scale (microscopic*) version of Maxwell 


equations to the macroscopic* (d -scale) equations of the electrodynamics of continuous media, 
taking averages in ascale a’ © Jad i 

This averaging method cannot lead to a bidomain model, no matter how the weight function is 
chosen for the averaging process in the a” -scale. In this sense bidomain models are not fully 


compatible with any a” -scale (microscopic*) averaging, as remarked by some authors [49]. 


The way out is to choose a different method of averaging that introduces three superposed 
continua: an intracellular phase, an extracellular phase, and a membrane continuum. The special 
procedure of phase averaging explained in Section (D), allows the construction of so-called 
averaged point fields, superimposed on the entire volume occupied by the tissue. From this 


construction it is possible to obtain the bidomain equations and fields as (,-scale phase averages 


of the already bulk d -scale equations and fields. These d -scale equations and fields are obtained 


from a’ -scale averaged equations and fields posed initially at the a -scale. 


(F) Choice of scales, averaged equations of continuity, membrane currents, averaged 
versions of Ohm law in cardiac tissue, bidomain boundaries and boundary 
conditions 


To begin, let's define the numerical values of the three scales assigned to cardiac tissue. 
From available histological data and from the results of the discussion done in section (D), we 
choose as a representative value of the microscale length d ~ 1m, both for the intracellular and 


for the extracellular spaces. 


22 The EQS approximation can be applied to study the fields originated in spontaneous physiological or 
pathophysiological processes, or to the production of electric currents in the tissues, in both cases using 
contact electrodes. The MQS approximation can be applied to mathematical models that describe the 
magnetic stimulation of biological tissues or the detection of the magnetic fields produced by these tissues. 
In both cases the magnetic measurement or the magnetic stimulation is done by non-contact devices. 
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As explained in section (C), the determination of space constants from electrotonic polarization 
of myocardium suggests choosing L ~1000 um for the characteristic length of the macroscale. 


Applying the geometric mean rule (equation (1) of section (D)) to these two scales, an estimation 
of a possible the mesoscopic length scale is ¢, ~ 33 um.” This will be the length of the sides of 


the cubic REVs used to take volume averages. 
We have dz=lum<«/, =33um<«L=1000um 
Equation (1) of section (D) express the approximate equality between the volume average and the 


area average of a smooth enough scalar point field y(t,y) unless a relative error 


y 2 
of . = o| mt This approximate equality can be applied to the cartesian components of 
smooth enough vector and tensor fields, and consequently to these two kinds of fields. 


So, the possibility to equate volume averages and the area averages over a representative 
p 2 L 
elementary volume (REV), unless a relative error 0} . | that is very small if - <1, justifies 


Ba 


the assumptions made in bidomain models about working with volume averages instead of area 
averages when dealing with electric current densities and the phase fractions associated with these 
currents. 

Once admitted the approximate equality of volume and area averages taken over a REV, it is 
straightforward to derive equations (5 a) of section (B): 


-V-J,=V-J,=4,J, 
Taking volume averages in equation (9) of section (E) and considering the electroneutrality at the 
mesoscale of both, the intracellular and the extracellular spaces, we obtain (V , j (t , x) =0 for 
the intracellular space and (Vv h(t )) =O for the extracellular one. *4 
Applying equations (23 a) and (23 b) to these last two equations we obtain, being 7,,, (x ) the 


membrane fraction in the REV, J,, (t ; x) the average membrane current there (positive if directed 


towards the extracellular space) and by definition J, = (i:) and J,= (i.) 


O=(V-7,\(t,2)=V- J; (t,2)+ %,(¥) J, (6.2) 


O=(V-7,)(t.2)=V-J, (t,%)- 2, (¥) Jn (t.¥) 


23 The mesoscale length must be larger than the linear dimension of the interaction region (including the 
damage region) of a microelectrode with the surrounding tissues. The estimation of this last length is usually 
less than 10 42m so everything is fine from this point of view. 


4 To obtain an emergent form of Ohm's law in each of the bidomain continuums, it is necessary that the 
transient phenomena associated with the charge transport processes relax towards electroneutrality on a 
time scale much smaller than the time scale characteristic of the variation in stimulating or spontaneous 
currents in the tissue. This is the same requirement needed to obtain the mesoscale equations of continuity 
in each continuum of the bidomain. 
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The meaning of J, (t ; x) is not as evident as it could seem to be at first sight. 


At the mesoscale it can be decomposed in a displacement and a conduction current like in equation 


(7) of section (B): J,,(®,,.W) = CA Oy Sin (©,,,W) 


However, to justify this last equation from considerations done at the microscale is much more 
complicated. 


Figure 14 shows a sketch of a hexagonal membrane unit, a concept apparently introduced by 
Kenneth Cole [4] [52]. 


Figure 14. Cole’s hexagonal membrane unit. On the average, only one ionic channel appears in a 
hexagonal unit. 


On the average, there is only one ionic channel per hexagonal membrane unit. 
So, the ionic conduction current must flow nearly parallel to the membrane surface, until an ion 
can enter to one channel that allows it to go through. 


related with the distributed membrane 


m 


The rest must be displacement currents Cam 


capacitance C, per unit area and with the velocity of variation of the local voltage D m defined 
at the mesoscale. 


The total ionic current i, (t 5X ) through the membranes of a REV with centroid at X can be 


ionic 
expressed as a sum in terms of the random currents [ x,, through the ionic channels that appear 


in these membranes: 


bionic x)= ¥[E4.] (1) 
k=1 \ [=1 


There are P types of ionic channels (perhaps a hundred different types) and there are ¢, channels 


of type k in the membrane surface A, (x) in the REV B (x) : 


If the REV is cubic and has a side length of ¢, ~ 33 44m, we can estimate an area of membranes 


S (A,,(%)) = 1000( um) . Considering the hexagonal unit, we can expect at least one ionic 


channel each 4x10~% ( um) of membrane. This gives 2.5x10° channels in the REV 


46 


August 2021 


membranes. If there are a hundred of different channel types, on the average there will be 
2.5x10° channels of the same type. 


If we assume that all the channels of type k are a sample of a statistical ensemble and if E [7 | 


is the mean value taken over the ensemble, then by the strong law of large numbers *°, almost 
dk P 

certainly we have: Mla See E[I, | and then we have I,,,,;. (t,x) & dif E{[I, | 
I=1 s 


Here q, is the number of channels of type k in the membranes of a REV. 
However, if 1, is the density of channels of type k (the number of channels of type per unit area 


of membranes in the REV), the g, =n,S (A ) so the density of ionic currents in the REV 


m 


membranes will be given by: 


bionic (t% Sen) (2) 


The expected values of the ionic channel currents are given by the kinetic equations of the 
membrane model chosen to describe the excitable membranes of the heart. These mathematical 
models range from the highly simplified model of Aliev and Panfilov to the highly detailed 
models that appear in the last editions of Cardiac Electrophysiology: from cells to bedside. 

In general, the kinetic equations of all the membrane channel types can be written in terms of a 


vector of channel gate variables W and a local membrane potential that must be interpreted as an 
area average [5] [8] [52] [62]. 

In phenomenological bidomain models, the transmembrane potential that appears in the channel 
kinetic equations is interpreted as the difference between the two electric potentials of the 


superposed continua, O, and O, as appears in equation (4) of section (B): 
®,,(t,¥)=®, (t,x) —®, (t,x) 
Let us analyse this mesoscale assumption working from the microscale. Analogously to was 


explained in section (D), for any point 5 €A,, cay the extension of the local point electric 
potential from the points of the intracellular (extracellular) phase ©. (Q .) of the REV B (x) ; 
toseA, (x ) in the phase boundary, is defined by the following limit: 


Pint S)= lim e(4F) © Pem(5) =, lim (5) ) 
YeQ; 0 B(x) FeQ, A B(X) 


An analogous definition applies to the cartesian components of a vector or a tensor field in both 
phases. 


The local point membrane potential 9, (t,5') is defined by 9,, (1,5) =9,,,,(t5)-Q,.,(t5) 
Its average (g,,), (t,x) over the membranes A, (x) in the REV B (x) is: 


°5 The probabilistic laws of large numbers are detailed, amongst other books, in Moran P (2002) An 
introduction to probability theory, Clarendon Press, Oxford, UK. 
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7 1 : 
Eaenelr Oe a, 5 


If the area averages of PD init ; 5) and @, mi Ss ) over A m3 ) are approximately equal to the 


t, 
volume averages of 9 (t.¥)=9(t, 9), (9) and 9, t,y)=9(t, ¥)1.(¥) over B(x), 


respectively, ae 
7 Lelenaa@r-7 { ele.s)aa(s=(oNe3)-(0NMoa) 
AS consequence be averaged membrane electric potential (Om) a (t,x) that appears in the 


kinetic equations of membrane models can be approximated by 
(Pm)a, (6%) *(Pida, (6¥)—(Pe)a, (6%) = 
D, (t,X)—®, (t,%) =©,, (t,%) 


The averaged density of ionic current through the membranes can be written now as a function 


(5) 


tOnIC 


eee (®,,,W) of a vector of channel gate variables W and a local volume averaged membrane 


potential D (¢ ae So, the mesoscale approximate decomposition of the membrane current 


density in a displacement and a conduction current density, like in equation (7) of section (B), is 
fully justified. 

From the mathematical tools developed in section (D) to be applied in volume averaging of fields, 
and from the quasi-stationary equations of electromagnetic theory for continuous media reviewed 
in section (E), the macroscopic versions of Ohm law for the intracellular and the extracellular 
continuum are easily obtained by volume averaging of the microscopic equations over the REV 
family. 

The volume averaged versions of equation (17) for the density of electric current in bulk are, with 


E(t,¥)=-Vo(t,¥): 


For the intracellular continuum: 


J,(t,8) =(Fouea)(-2)=-(6;,-Ve,)(t,8) +(6;,-Ep, }(t,2) (6a) 
For the extracellular continuum: 
F,(t,8) = (Grate (0%) = (6, °V@,)(t,8) + (6, Ep. )(t,%) (6b) 


Applying to these continua the results obtained in section (D) for the averaged transport of a 
substance-like quantity driven by a local point scalar potential, we obtain: 


(&,-VQ,)(t,¥) = f, (¥)G, (X)-V®, (t,x) (7a) 
(6.-Vo,)(t8)= f.(8)6,(3)-VO, (13) (1) 
Here, by definition, ,(t,%)=(g,)' (t,¥) and ®,(t,x)=(g,)! (t,x) are intrinsic phase 
averages, taken over the REV with centroid X, of the local point electric potentials 
9, (t, ¥)=@(t, ¥)1,(¥) in the intracellular space (phase i) and @, (t,¥)=9(t, ¥)I, (5) in the 
extracellular (interstitial) space (phase e ). As always, I, (¥) and I, (y) are the phase indicator 


functions defined by equation (4) in section (D). 
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The effective macroscale conductivity tensors G, (x ) ‘ G, (x) are derived from the microscale 


conductivity tensors G, (¥) and 6, (¥ ) by one of the two closure methods developed in section 


(D). They are symmetric and definite positive. 


In one of this closure methods, a reciprocity theorem and a closure hypothesis are applied as 
shown in section (D) for the uncoupled linear transport of a substance-like quantity, like electric 
charge. 


In the other, when applied to cardiac tissue, a representative unit element (RUE) of the tissue’s 
microstructure is introduced. 

In general, this RUE does not coincide with a REV. It is a kind of statistically average of the 
REMs defined in a relatively homogeneous region considered from the macroscale. 

When the tissue shows non-homogeneities at the macroscale level, a minimal set of suitable 
defined and representative RUEs must be introduced to approximately describe the macroscale 
spatial variation. 

Assumption A made in section (D) about the difference of flow resistance between the bulk of a 
phase and its interface, allows to impose as a first approximation, zero normal flow through the 
interface. 

This same assumption is also justified when the interface corresponds to the biological 
membranes separating the intracellular from the extracellular space. 

Let us consider an experiment with controlled voltage [Cole]. 

Suppose that a charge equivalent to 200 monovalent ions cross the hexagonal unit of membrane 
each millisecond. 

During this same interval of time, each side of this hexagonal unit suffers the impact of 6x10 ° 
ions. 

So, only one ion each 3x10 7 can cross the hexagonal unit: this is a formidable resistance to current 
flow. 

With this remark, we see that the same procedure developed in section (D) can be implemented 
in the case of myocardium to derive the effective conductivity tensors, solving the closure 
problem described there, in each selected RUE. 

The effective (macroscale) conductivity tensors are defined by the following formula (formula 
(54) of section (D)): 


Fy (8) G2) = (4 (5)-(7+ VB, (3)))() u=ie 
The vector fields b, ( y ) and b, ( y ) are solutions of the elliptic boundary value problem given 
by equations (48), (49), (50) and (51) of section (D), with K, (¥)=6,,(¥) and the Gray 
decomposition G,, (¥) = (6, \! (x) + o, (Cy) : 
In the corresponding phases ( 4 =i,e ) of the RUV the fields b i (¥) verify the equation: 
V-(6,,(5)-Vb, (5) =-V-(4,/ (5) =i,e 
To this, we must add the boundary conditions in the membrane interface, the periodic conditions 


in the RUV boundary and (b, ) (X) =0 for both fields. 


In electrophysiology, generally the Planck field can be neglected when the electric current in bulk 
is described. 

If this assumption is also done here, from equations (6) and (7) we derive the effective Ohm laws 
for the intracellular and the extracellular continua that appear in section (B): 
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1,(t.8)=-f(8)6,(8)-VO,(48)— F.(48)=£,(8)6.(8)-VO, (1.8) 


Now, let us consider an operational procedure suitable to define a sharp bidomain boundary (at 
the mesoscale) from the fuzzy distribution of structural details and physical properties at the 
histological level (microscale). 

Histologically, there is no real solution of continuity between the extracellular space and the extra 
tissular space. 

In addition, the membranes of the outer myocardial cells present numerous microscopic details 
(invaginations, folds, and other structural details) that are not relevant at the level of averaged 
fields. By defining the boundary of the biological tissue at the mesoscopic level, these details 
must be eliminated. 

Suppose that we start from a point P located inside the myocardium and advance on a ray 
emanating from this point, as shown in Figure 15. 


Figure 15. A beam of rays emanating from a point within the biological tissue. 

Along this ray we move the centroid of a REV and for each position we find the volume fraction 
of the intracellular space. 

If we introduce an abscissa s originating at point P, the volume fraction f, (s) of the 


intracellular space varies between certain characteristic values of the tissue interior (depending 
on the anatomical characteristics and the physiological state) and zero, when the centroid of the 
REV is far enough from the heart tissue, inside the exterior region signalled with o in figure 15. 


Figure 16 shows the behaviour of /f, (s) as _s increases. 


{6) 


Ss 


Figure 16. The volume fraction of the intracellular space as a function of the distance along a ray directed 
towards the extra tissular space. 
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On a segment of values of the abscissa, intermediate between the interior of the myocardial 


tissue and the extracellular space, the function f, (s) varies rapidly, on a scale comparable to 


L 


the mesoscopic scale £,, or less, instead of varying on the macroscopic scale L as occurs well 
in inside the myocardium. Inside this interval of fast variation, the boundary of the tissue must 
be located. 


One possible method of constructing the tissue boundary is to select an auxiliary surface within 
the tissue, as suggested in Figure 17. 


Figure 17 The auxiliary surface and the beam of rays constructed from each point of the auxiliary surface. 


From each point several rays are drawn directed towards the outside of the biological tissue. 


Figure 17 shows this construction for three points 7, , P, and P, . From each point, the ray that 
presents the sharpest variation of ff, (s) with distance is selected. Then the point of that ray is 


determined at which f, (s) takes half of its average value in the last segment prior to the abrupt 


variation segment (see Figure 16). The points thus determined define the boundary of the 
biological tissue. 

It is convenient to choose the auxiliary surface relatively close to the "blurred area" which, from 
the histological point of view, constitutes the transition zone between the interior of the 


myocardium and the extra myocardial space. In the "blurred zone" f, (s) varies, in a normal 


direction to the boundary thus constructed, so rapidly that it can originate a spatial scale smaller 
than the mesoscale. Here the relationship of scales on which the theory of averaged fields rests is 
lost. We have a true "geometric singularity" within the framework of this theory. This problem is 
solved, as explained above, by introducing an ideal, regular boundary surface, where the real 
boundary is, in a sense, indefinable and where the micro geometric details suggest a description 
by an opposite idealization: a fractal surface. 


The value of a volume averaged field in each point of the idealized boundary will be defined as 
its value when the centroid of the REV is located at this point. 


The above suggested method for boundaries definition at the mesoscale level can be applied also 
inside the heart tissue to separate infarcted and fibrous areas of the heart muscle from healthy 
myocardial tissue. 
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Once constructed a mesoscale boundary between the cardiac tissue and an external continuum, 
and under the same assumptions made in section (B), then equations (16) and (17) can be applied 
to describe mesoscale boundary conditions in the bidomain model. 


(G) Two continua: passive and active bioelectric continua 


In this section the bidomain equations are rewritten to introduce two new continua, both 
heterogeneous and anisotropic: a passive bioelectric continuum and an active bioelectric 
continuum. These media are coupled by two mechanisms: a surface coupling at the bidomain 
boundary and a volume coupling through a term that generalizes the activating function to 
anisotropic syncytia. 


Let us begin with the set of bidomain equations for cardiac tissue. 


The volume averaged continuity equations. These relate the local volume averaged density of 


electric currents in the intracellular continuum J; (t es ) and in the extracellular (interstitial) 
continuum J, (t ; x) with the local averaged membrane electric current density J,,, (t Xk ) and the 


local averaged membrane fraction 7, (x) : 


V-J,(t,X)+ x,, (%) J, (t,¥%) =0 (la) 
V-SI.(t,¥)— Xm (X) J, (t,%) =O (1b) 
The macroscale Ohm laws relating the local volume averaged current densities with the gradients 


of the intrinsic local volume averaged electric potentials, the volume fractions of each phase and 
the corresponding local macroscale conductivity tensors are: 


J, (t,%) =—f, (%)G, (X)-V®, (t,x) (2a) 
J ,(t,x)=—f, (%)G, (¥)-V®, (t,x) (2b) 
The volume fractions of each continuum are the integrals of the indicator functions over the REV: 
f(S=s J u(s)av(3) Ga 
B(x) 
£(@== J LG) av(5) (3b) 
B(x) 


The macroscale intrinsic local phase averages of the microscale electric potential are: 


e(ni)=z J oC NLE)AV() (4a) 
©. (n8)=— J o(t.3)1,(5) av (3) (4b) 
e B(x) 


The local volume averaged membrane voltage given by: 
®,, (t,X) =, (t,%)-®, (t,x) (5) 
The equations that describe the local averaged membrane electric current density (when the 


electric consequences of membrane invaginations are neglected) and the kinetics of its channel 
variables: 
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J,( pW) =C, 20 4D ge (Oyo) (6x) 
<W = F(,,W,3) (6 b) 


These equations are verified for every time instant and in every point of the bidomain {that 
represents the electric behaviour of cardiac tissue at the macroscale level.”° 


To apply these bidomain equations in a theoretical approach to the excitation and propagation of 
an action potential in cardiac tissue, due to the stimulation produced by pacemaker electrodes, it 
could be convenient to substitute the two superposed intracellular and extracellular (interstitial) 
continua by a new continuum, the passive bioelectric continuum, superposed with the 
continuum of excitable membranes and electrically coupled with this last one. When considered 
in conjunction with the passive bioelectric continuum, we call the membrane continuum active 
bioelectric continuum. 


The local volume averaged electric current density J (t, x ) that flows in cardiac tissue is the sum 


of the current densities (2a) and (2b) in the intracellular and extracellular continuums: 


71,8) =F, (2)G, (2)-V®, (1.2)-F,(2)G.()-VO,(42) 
An inspection of equation (7) suggests that in general there is no possibility to derive from it an 
emergent Ohm law to describe the global electric current flow in the mesoscale level using the 
gradient of an average electric potential: 


O(a J (6.5) AV (= HH B(3}+ £(5)&.(63) (8) 


Only if the volume fractions are constant and the conductivity tensors of the intracellular and 
extracellular continuums are equal, an Ohm's law results in terms of the gradient of potential (8). 
But it is trivial and has no interest in the case of electrical current flow in the myocardium. 


However, an auxiliary scalar field Y (tj x) can be constructed occupying the whole bidomain 


region. 
Define first a global, symmetric, and definite positive (and consequently invertible) conductivity 
tensor: 


G(¥)= f,(¥)G,(¥)+ f.(2)G. (%) (9) 
G (x) corresponds to two contributions in parallel to the global current flow in the biological 
tissue, following a kind of mixture law applied to the results of an averaged field theory. 7” 
Then the auxiliary field W (t,%); called the electric pseudopotential is introduced by the 


following equation: 


VW (1,¥)=G"(x)-(f,(¥)G, (¥)-V®, (.¥)+ f, (XG, (¥)-VO, (1,3) (10) 


6 The explicit dependence of the position vector in the bidomain allows to describe at the macroscale a 
certain degree of membrane heterogeneity. 


7 The relation of the different versions of the law of mixtures and mean field theories can be found in the 
book by Snarskii A, Bezsudnov I, Sevryukov V et al (2016) Transport Processes in Macroscopically 
Disordered Media: from Mean Field Theory to Percolation, Springer, NY USA. 
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From equation (7) and (11): 

J (t,¥) =-G(X)-V (t,x) (11) 
From equations (1) and (5) we have V - J (i; x) = (0 in the bidomain. 
From this last equation and from (6) we obtain that the electric pseudopotential verifies the 


following Laplace-like elliptic equation in the region $R occupied by the bidomain: 
V-(G(#)-V¥ (7,%))=0 (12) 
If the normal component of the global current density 7- J (4, x) = h(t, x) is known in each 


point X the boundary O® and for every instant of time, then the following non homogeneous 


Neumann boundary conditions are added to equation (11): 
i-(G(%)-V¥(t,¥))=A(t¥) V¥EOR (13) 


Equations (11) and (12) describes the passive bioelectric continuum. 


Now, let us describe the active bioelectric continuum. 
From (5), (10), defining A(x) =G(X) ' - f,(¥)G,() and B(x) =G(x) '- f.(¥)G, (%) 
we obtain: 
VO, (t,x) =V¥ (t,x)-A(X)-VO,, (¢,%) (14a) 
V®, (t,%) = VY (t,x) + B(X)-V®, (t,x) (14b) 
From equations (2a) and (2b) and defining C(x)= if (X)G, (x): G(x)" -f, (#)G, (x) it 
follows: 
—J,(t,%) = f (¥)G, (%)-V P(t, ¥)+C(X)-VO, (t,x) (15a) 
—J,(t,%) = f, (¥)G, (%)-V¥(t,%)-C(Z)-VO, (¢,%) (15b) 
Adding member to member equations —V-J,(t,x)=7,,(X)J,,(t,%) and 
Ved, (Ba) Ln (x) Tin (ia), considering equations (15), we obtain this equation: 


wear a 1 te ia Hes wat ech = i >. 

V-(E(3)-70, (3) + 3 £(3)6(3)- £.(3)6,(3))-V¥(43))= Nee) 
From this last equation jointly with equations (6a) and (11), and defining 

£ 1 Z 4 , 

D(%)=5-(£.(2)G.(8)-£ (2)G, (2) )-E"(%) (16) 
we have: 
V (C(%)-V®,,)+V (D(x): (t,x)) 29. (2 Ne <9, tT onic (0..")] (17) 
This equation jointly with equation (6b) that gives the kinetics of the channel variables, describes 
completely the active bioelectric continuum. 


The term V- (D (x ) J (t oe )) is a generalized activation function for the bidomain models when 


considered from the point of view that decomposes the continua in a passive and an active 
bioelectric one. 


It couples the active bioelectric continuum (characterized by the state variables D,, ,W ) with the 


passive bioelectric continuum, in this case represented by the global density J (t, x) of electric 
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current flow in the biological tissue. This last field can be obtained independently, considering 
equation (11) and solving the elliptic boundary value problem given by equations (12) and (13), 
that describe the passive bioelectric continuum. 


Now, consider the meaning of the pseudo potential introduced to describe the passive bioelectric 
continuum. 


Let us suppose that the cardiac tissue is at rest until the time instant t=O. 
As consequence the membrane voltages have their rest values everywhere. 


From t=0 upto t=t » arectangular pulse of controlled current is applied by external electrodes. 


Due to the Ohmic behaviour of the passive bioelectric continuum, immediately after the pulse is 
applied the pseudo potential undergoes a step variation in each point of the bidomain, as suggested 
in Figure 18. 

Of course, these variations are position and pulse amplitude dependent. 

As the point averaged membrane potential cannot undergo a sudden variation, the point averaged 
potentials in the intracellular and extracellular spaces, at each point of the bidomain, must undergo 
the same sudden variation, at the beginning of the pulse. 

As long as the external pulse of constant amplitude lasts, the pseudo potential at each point must 
remain constant. 

However, the membrane potential leaves its rest state while the membranes are polarized. 


Fi, Be Yr 


Figure 18. The time variation of the pseudo potential ¥ (t ix ) and the point averaged potentials D j (t 8 ) 


in the intracellular and , (t JX ) in the extracellular continuum, in a fixed centroid X ofa given REV 


B(x): 


Due to formula (5) and the dynamics of charge transport in the intracellular and intracellular 
continuum, this variation in the membrane potential must imply a time variation of both, the 
intracellular and the extracellular potential, as is suggested in Figure 18. 

When the external current pulse ends, the pseudo potential drops to zero everywhere, as does the 


global density current J (t, x ) on a macroscale, averaged over the REVs. 


But this does not mean that the point currents densities averaged over REVs in the intracellular 
and extracellular continua disappear, but rather they adjust the values that they adopt so that the 
macroscale global current density remains equal to zero. 

At the end of the pulse, the intracellular and extracellular point potentials averaged over the REVs 
undergo the same jump at each point of the bidomain, but in the opposite direction, as the jumps 
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they had suffered at the beginning of the pulse. After the end of the pulse, they evolve tending 
asymptotically towards zero, as suggested in Figure 18. 

While the intracellular and extracellular potentials of the bidomain are affected by charge transfer 
across cell membranes, the pseudo potential is not directly affected by these processes. 

In accordance with these considerations, the value adopted by the pseudopotential when a 
rectangular pulse of controlled current is applied by external electrodes, coincides with the 
variation suffered by the point averaged potentials in the intracellular and extracellular continua 
when the pulse begins. 

But although the pseudo potential is not a true electric potential, it is possible to think of it as the 
electric potential of an equivalent Ohmic volume conductor, inhomogeneous and anisotropic, 


whose properties are characterized by the tensor G ( x ) = i (x ) G. (x ) ad i F (x ) G, (4 ) : 


To end this section, let us see possible applications of this alternative approach to bidomain model 
in electro cardiology. 

In references [53], [54] and [55] the above derived coupled equations for the passive and active 
bioelectric continuum were applied to a physical-mathematical analysis of the excitation of 
myocardium by pacing electrodes with focus on the just-threshold dynamics of the electrode- 
myocardium system. 

A purely analytical approach was developed, based on nonlinear modal analysis techniques, 
without resorting to ab-initio methods of digital simulation of dynamics. 

Several experimental results are thus explained by the obtained analytical formulas. 

Among the results of the theory, are the apparition, size and shape of different liminal regions, 
and the behaviour of the parameters of strength-duration curves when the size, shape and distance 
of the pacing electrode, or the properties of biological tissues are modified. *® 

Some new facts are predicted, mainly related with liminal regions, the family of electrode- 
myocardial threshold states and the parameters of the strength-duration curves. 

This theoretical approach can be extended to other excitable tissues, such as electrically coupled 
smooth muscles. 

Besides, it can be (and was) used as a guide to the design of pacing electrodes. 


(H) Discussion and conclusions 


A solution was given here to the four main problems that were distinguished in relation to 
bidomain models at the end of section (B), including the error estimations of the derived formulas. 
The introduction of the three superposed continua, including the fields that appear in the equations 
of bidomain models, was justified considering the complexities of tissue’s microstructure at the 
histological level. 

The continuity equations and the macroscopic Ohm's laws that are postulated in the bidomain 
theory were obtained here from the equations that govern the local point fields on a microscopic 
scale, using tools from the theory of averaged fields. Parameters such as the fraction of membrane 
area per unit volume or the conductivity tensor assigned to each domain, result from a 
combination process, at the mesoscale level, of irregularities in space and time, which occur at 
microscopic scale associated with tissue microgeometry. 


8 Suarez-Antola, R, The time constants for cathodic make stimulation of electrical syncytia: an application 
to cardiac pacing, Proceedings of the 2006 Conference of the IEEE Engineering in Medicine and Biology 
Society, 4031-4034, 2006. 
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The definition of the transmembrane current density as a local averaged field makes it possible to 
connect the random currents in the ion channels with the kinetic models of the membranes of the 
different types of cardiac cells, considering that these models are averages at the mesoscopic level. 


The averaged local transmembrane voltage appears as an independent concept from the definition 
assigned to it in the bidomain theory, but the equality between the two is demonstrated in the 
framework of the theory of point averaged fields, from the relationship between averages over 
areas and averages over volumes. 


The border between the bidomain and the extra-tissue space can be defined by taking as a guide 
the volume fraction corresponding to the intracellular continuum. It corresponds to a region in 
which the relationship between the micro, meso and macro scales is destroyed. In this area the 
spatial scale of variation of the field of intracellular volume fraction changes abruptly from being 
a magnitude of macroscale order to being a magnitude of mesoscale or even microscale order. 
This constitutes a singularity from the point of view of the averaged fields and justifies the 
introduction of a boundary surface of the form that was described in this work. 


Even when the conductivity properties at the microscale level were those of a homogeneous and 
isotropic conductor, the geometric impediment to the flow of electric current in the interstitial 
continuum or in the intracellular continuum due to the irregular distribution of the membranes on 
a microscopic scale can be expressed, on a mesoscopic scale, by means of a tortuosity vector. 
This vector is analogous to that used to describe the geometric obstruction to mass-diffusive 
transport in the pore space of a saturated porous medium. 


Two new continua were proposed to reformulate the equations of bidomain models in electro 
cardiology: the passive and the active bioelectric continuum. 


The generalized Laplace equation derived for the pseudo-potential corresponding to the passive 
bioelectric continuum can be solved with suitable boundary conditions totally independent of the 
field variables corresponding to the active bioelectric continuum. 

In the study of myocardial excitability by an extra tissular active electrode, it is convenient to start 
by finding the pseudo-potential and the current density field generated by the electrode in the 
passive bioelectric medium. 

Once the global electric current density field has been found in the passive bioelectric continuum, 
it is possible to face the resolution of the field equations for the variables that characterize the 
active bioelectric continuum. In these equations a term appears that depends on the global current 
in the passive bioelectric continuum, through which the processes in the active continuum are 
coupled with the processes in the passive continuum. 


The uneven anisotropy characteristic of myocardial tissue results in the fact that, even when the 
conductivity tensors do not vary with position, the current density in the passive bioelectric 
continuum is coupled at all points with the processes in the active continuum. This fact makes the 
depolarization and hyperpolarization of the active continuum membranes present different 
characteristics from those observed when the anisotropy is equal. 

In the framework of the original equations of the bidomain models of electro-cardiology, these 
motivated numerous papers related with cathodic make, anodic break, cathodic break, and anodic 
make excitation of myocardium, including a comparative analysis of similitudes and differences 
with nervous tissue [8] [62]. 


To clarify the relation between the bidomain models of electro cardiology and the classical 
electromagnetic theory of continuous media, it was necessary to distinguish two very different 
levels of averaging. 
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One level of averaging allows to pass from the a -scale (microscopic*, a = 0O,1nm) version of 


Maxwell equations to the macroscopic* (d -scale, 142m) equations of the electrodynamics of 


continuous media, taking averages in a scale a’ © Jad =10 nm. 

This averaging method cannot lead to a bidomain model, no matter how the weight function is 
chosen for the averaging process in the a” -scale. In this sense bidomain models are not fully 
compatible with any a° -scale (microscopic*) averaging, because they cannot be obtained by 
simply increasing the size of the volume used for averaging Maxwell-Lorentz equations. 


We found that the way out is to choose a different method of averaging that introduces three 
superposed continua: an intracellular phase, an extracellular phase, and a membrane continuum. 


A family of restricted representative elementary volumes of linear size /,¥./L.d ~33 my is used 


to construct the volume averages. The special procedure of phase averaging allows the 
construction of so-called averaged point fields, superimposed on the entire volume occupied by 
the tissue. From this construction it is possible to obtain the bidomain equations and the 


corresponding bidomain fields as (,-scale phase averages of the already bulk averaged equations 


and fields, obtained from a” -scale averaged equations and fields posed initially at the a -scale. 


The scale L of significant spatial variation of the point averaged fields is of the numerical order 
of mm. 


Finally, the following remark. Let us assume the new explanatory framework based in non- 
equilibrium phase transitions in the cell membrane and its adjacencies [40] [41] [42] could gain 
importance with respect to the current explanatory framework based on a multiplicity of 
molecular receptors, channels, and pumps. Then, neither the main aspects of the bidomain models 
in electro cardiology nor the type of foundation of these models given here should be significantly 
modified (apart from applying the law of large numbers to the ensemble of channels in the 
membranes of a REV). 


(IT) Appendix: Phasor vector field version of Maxwell equations corresponding to 
monochromatic oscillations, and a generalized perturbation approach suitable to 
describe biological tissues. 


Let us consider again the version, due to Minkowski, of the Maxwell’s equations for a continuous 
medium: 


~ . @D . 
VAH =j+— 1 VAE=-— p 
\ J Ot (1) x ar (2) 


V-B=0 (3) V-D=p (4) 


Here, as usual, H (& y) represents the magnetic field, B (t, r ) the magnetic induction field, 
E (@ y) the electric field, D (t, y) the displacement field, J (t, y) the electric current density of 


free charges, and 9 (t, iz ) the density of free charges. 


Taking the divergence in both members of (1) and considering (4) the following version of the 
conservation of electric charge is derived: 


V-j+—=0 (5) 
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It is possible to study the fields near the stationary limit, remaining in the time domain. To do 
this, the time variable is redefined, introducing a frequency as a new independent variable, and 
developing all fields in a power series of this frequency [23]. 7’ 


To order zero, we obtain the uncoupled equations of the electrostatic and magnetostatic fields: 
VAHO=79 VAEO=60 V-B!=0 VD) v-j=0 6 


The fields j and p”) are assumed to be known. 


In general, the approximation to order m is: 


z 7 > (m-1) . (m1) 

VA” =704 2 yng BBM 0 VD = po” 
Ot Ot 

V-j cr 226 MAIS: (7) 


Ot 
The fields of order m are related with the fields of order m-1 in such a way that fields of all orders 
can be obtained in succession from zero order fields, adding suitable boundary conditions (that 
are developed also in powers of the frequency) and posing initial values [23]. 


The quasi-stationary approximation in the broad sense consists of zero order and first order terms 
only. The quasi-stationary approximation in the strict sense consists of zero order terms only. 

In the broad sense, it gives a good enough description of relatively slow phenomena such as those 
observed in many electric circuits and in the electric or magnetic stimulation or recording in 
biological tissues. 

As the speed of the phenomenon of interest increases, the higher order terms in the approximation 
becomes, progressively and in succession, more and more important. 


How do we measure slow or fast? A suitable measure is to divide the transit time /,,, that an 


electromagnetic perturbation needs to traverse a characteristic length l . of the studied region (for 
example, the size of the volume conductor formed by the cardiac tissue) by a characteristic time 


of variation of the signal of interest f.,, (for example, a propagated action potential in the heart). 


l 


Given that c is a representative speed of the perturbation, then ‘,,, = — and dimensionless 
6 


t 
number €= — can be used as a measure of the degree of speed of an electromagnetic 
ch 
phenomenon in relation to the geometric and physical characteristics of the system in which it 
occurs. The smaller is this number, the slower is the phenomenon. 
When the medium is linear, nondispersive, non-homogeneous and isotropic, it is straightforward 


to rewrite Maxwell equations in dimensionless form. *° 


t 
In the Faraday law and in the Ampere-Maxwell law the small parameter ¢=—~ appears 
ch 
multiplying the terms with partial derivatives relative to dimensionless time. A development of 
all the dimensionless fields in powers of the small parameter gives a sequence of successive 
approximations equivalent to equations (6) and (7) given above. 


?° This approach was developed by Fano, Chu and Adler in their book Fano R, Chu L and Adler R (1968) 
Electromagnetic fields, energy and forces, MIT Press, Cambridge, Mass USA. 

3° The details of the construction and a detailed analysis of the dimensionless equations for a linear, non- 
dispersive and isotropic medium can be found in [18]. 
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However, the cardiac tissue is anisotropic and highly dispersive. As shows Figure 13, the global 
relative dielectric permittivity varies more than two orders of magnitude between the band of 
frequencies of interest here, comprised in the interval 0 tol0 kHz. 

Due to this, it seems better to work with another approach to study all the fields near the stationary 
limit: an analysis in the frequency domain. 


So, let us consider the vector phasor field solutions of Maxwell’s equations corresponding to 
monochromatic oscillations of angular frequency @ [43] [47]. 


We substitute the fields E(t,¥), D(t, ¥) : H(t,¥) and B(t,y) by eE(o:3), 


e!"D(a; y) : eH (a; y) and e!"B(a; y) respectively, the conduction current I(t, y) 
by e!"J(@;y), and p(t, y) by e! (ay). 

Here E(a;y), D(a;7) ,H(a;y), B(a;y), J(@;y¥), and p(@;y¥) are phasor 
vectorial amplitudes (this is represented by the point that appears over each arrow) and, as usual 
in electric engineering, j =V—1. 

Substituting these vector phasor field solutions in Maxwell equations, we obtain the following 
equations in the frequency domain: 


> 


VAH(a;3)=J(@;3)+ joD(a;5) (8) 
V A E(@;3)=—joB(a;5) (9) 
V-B(a;%)=0 (10) 
V- D(a: 3) = p(a:5) (11) 
And the phasor version of the free charge balance equation: 
V-J (a3) + jop(a:¥)=0 (12) 


The following constitutive relations describe a linear, time invariant, anisotropic, non- 
homogeneous and lossy dispersive medium in the frequency domain: 


J (@:3) =F, (a; ¥)+6(@:%)-E (0:5) (13) 
D(a; 3) = &(@;3)-E(@; 3) + B (a; ¥) (14) 
B(o;¥)= 4H (a;¥) (15) 


The components of the dielectric permittivity tensor é (a; y) and the electric conductivity tensor 
6(a; y) (all of them defined in the frequency domain) are in general complex scalar fields whose 
components are the Fourier transforms of the component kernels €; (t, y) =O. ea xv ij (t ,y ) and 
On (t , y) of the causal integral operators [ E ] (t,¥) (equation (14) of section (E)) and 
6[E ] (t, ¥) (equation (12) of section (E)) respectively. * 


31 A presentation of the foundations of electromagnetism, including the analysis in the frequency domain, 
can be found in the books [20], [23], [47] and in the concise book by Post E (1962) Formal structure of 
electromagnetism, North-Holland, Amsterdam, Netherlands. 
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From now on, we suppose that they are symmetric. 

The term J, (a; ¥) is the Fourier transform of the electric current density ve its y) that 
encompasses the bulk effect of the concentration gradients of mobile ions in the fluid spaces 
delimited by the cytoskeleton and the micro-trabecular network, the entrainment due to the 
possible bulk flow of tissue fluid and localized ionic flows through open channels in cell 


membranes, driven by electrochemical potential differences and possibly related with 
nonequilibrium phase transitions in biological materials, as was said in section (E). 


The term P. (a; y)i is the Fourier transform of the polarization term P. (t, y) that encompasses, 


both the bulk effects and localized membrane effects, on bounded charges, related with 
conformational changes and other cooperative effects in macromolecular structures, as mentioned 
in section (E). 

According to (15), the biological tissue is considered as transparent to the magnetic field, in the 
sense that B(t, y) = LH (t, y). 


If é(a) and o(a) are, respectively, representative (maximum) values of the dielectric 


permittivity and the electric conductivity at each frequency: 


é (a; 3) =e(@)é (a5) (16) 
6(@;¥)=0(a)o(a;5) (17) 
The elements of the dimensionless tensors é (@; y ) and 6(a; y) are of the numerical order of 


unity or less. 


If ia is a characteristic length of the studied region, we introduce the dimensionless space 


variables Y =—— Y . Then, introducing the operator V defined by partial derivatives relative to 


Cc 


dimensionless spatial variables: 


V=—V (18) 


Next, we select a representative value F, of the electric field. In fact, this quantity is going to 


disappear during the process of obtaining the dimensionless equations. Then we put: 


E(o;3)=E,E(0:3) (20) 
We introduce a dimensionless displacement field: 
D(o;7)=«e(@)E, D(a;7) (21) 


Taking as reference the relation between the magnetic and the electric fields in the far field of a 
monochromatic radiating dipolar system in linear, homogeneous, isotropic, and non-conductive 
medium, we introduce the dimensionless magnetic and induction fields: 


H(o;3)=E,,|" ae 1) Fw y) (22) 


peer. ) Mo H (o; y) (23) 


We introduce a dimensionless electric current density and electric charge density: 


J (a; ¥)=o(@)E,J(@;¥) (24) 
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* 


aie ls c ee 
AO) ayn) (25) 


Finally, we define the following characteristic times and dimensionless numbers: 


Monochromatic propagation time: box (a) =Jé (a) Le a (26) 
Magnetic diffusion relaxation time: tin (a) = Ho o(a) 4 i (27) 
Electric charge relaxation time: t, (a) =— (28) 


A dimensionless measure of the degree of speed of a monochromatic electromagnetic oscillation 
in relation to the geometric and physical characteristics of the system, in which it occurs, at this 
same frequency: 


E om(@) = 1, (@) (29) 
The smaller is this number, the slower is the phenomenon. 


The dimensionless alfa number, and the relation between the three dimensionless times given by 
equations (26), (27) and (28): 


a(a) = nl) _ tol) | foo | le) l, (30) 
t (a) ein (a) & é. (a) 
Considering all the definitions introduced above, after a symbolic calculation process applied to 


the phasor equations (8) to (12), and equation (15) (magnetic transparency) we obtain the 
dimensionless form of Maxwell-Minkowski equations: 


=> 


VAH (0:3) =a(o)J (0:5) + 7 €,,(@) D(@5) G1) 
VAE(o;¥)=—j€.,, (@) H(a;¥) (32) 
V-H(w.5)=0 (33) 
V-D(o;%) = p(o: y) (34) 
VJ (0:9) +f Eom (a) p(@;,¥)=0 (35) 


The dimensionless form of the constitutive equations, applicable to cardiac tissue, are: 
J(a;¥)=J.(@;9)+6(@;¥)-E(@;¥) (36) 


D(a; 3) = P, (a; 3) + E(a;¥)-E(a;5) (37) 


To order zero in €,,, (a) we obtain this version of the strict quasi-stationary equations: 
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VAE «(a 3)=6 (38) 
V-(&(05)-Eo(@3)]= Bio 5)-V-B. (0:3) 39) 
VA Hy) (o:5)= a(0)(F.(o: 3)+6(a5)-E (a 5) (40) 


V-H (a5) =0 a) 


Here P,(@;y) and J (@; ¥) are considered as given. They appear only in the 


approximation of order zero. 


So, to the first order in €,,,, (a) : 


V AE, (0: ¥)=—j H)(@¥) (42) 

VAH, (0:3) =a(0)E(0:3)-E (09) + JDio(@:3) (43) 
V-(&(o3)-Eu(o: 5) = P(@7) (44) 
V-H,)(@;%)=0 (45) 


Adding the approximations of order zero and one we obtain the quasi-stationary approximation 
in the broad sense. 


€ em ( @) 


All this presupposes that the quotient of dimensionless numbers is small relative to 1. 


From equations (26), (29) and (30): 


= = ot, (a) (46) 


According to Plonsey and Heppner [35], the dimensionless number ot, (a) is less than 0.15 in 


the interval of frequencies comprised between 0 and 1 kHz. 
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